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The  research  has  been  carried  out  simultaneously  on  three  aspects  of  aircraft  wings 
performance  optimization,  as  follows: 

i.  Structural  tailoring  and  modern  control  of  thin-walled  model  of  composite  aircraft  wings. 

ii.  Structural  tailoring  of  a.  low-aspect  ratio  plate  model  of  composite  aircraft  wings, 
ni.  Integrated  structural  design  and  vibration  control  by  multiobjective  optimization. 

Significant  progress  has  been  made  on  all  fronts. 
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Abstract 


The  research  has  been  carried  out  simultaneously  on  three  aspects  of  aircraft  wings 
performance  optimization,  as  follows: 

i.  Structural  tailoring  and  modern  control  of  thin-walled  model  of  composite  aircraft  wings. 

ii.  Structural  tailoring  of  a.  low-aspect  ratio  plate  model  of  composite  aircraft  wings. 

iii.  Integrated  structural  design  and  vibration  control  by  multiobjective  optimization. 
Significant  progress  has  been  made  on  all  fronts. 
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Introduction 


The  research  can  be  divided  into  the  following  areas: 

i.  Structural  tailoring  and  modern  control  of  thin-walled  model  of  composite  aircraft  wings. 

ii.  Structural  tailoring  of  a  low-aspect  ratio  plate  model  of  composite  aircraft  wings. 

iii.  Integrated  structural  design  and  vibration  control  by  multiobjective  optimization. 

Five  technical  papers  have  resulted  from  the  research.  The  research  is  reported  according 

to  these  papers. 

1.  Refined  Structural  Modeling  For  Enhancing  Vibrational  and  Aeroelastic  Characteristics 
of  Composite  Aircraft,  Wings,  L.  Librescu,  L.  Meirovitch  and  0.  Song* 

An  analytical  study,  including  tailoring,  of  the  vibrational  and  static  aeroelastic  response 
characteristics  of  anisotropic  composite  aircraft  wings  in  the  form  of  thin- walled  beams  is  pre¬ 
sented.  The  theoretical  analysis  and  numerical  results  encompass  effects  such  as  transverse 
shear  flexibility  exhibited  by  the  advanced  composite  materials,  warping  restraint  charac¬ 
terizing  cantilevered  structures,  elastic  anisotropy  and  induced  structural  couplings.  The 
complex  effects  of  these  factors  are  highlighted  and  the  power  of  the  tailoring  technique 
toward  enhancing  the  dynamic  and  static  structural  characteristics  is  demonstrated. 

The  paper  has  been  presented  at  the  AIAA/ASME/ASCE/AHS/ASC  34th  Structures, 
Structural  Dynamics,  and  Materials  Conference,  La  Jolla,  CA,  April  19-22,  1993,  and  has 
been  accepted  for  publication  in  the  journal  La  Recherche  Aerospatiale.  A  copy  of  the  paper 
is  enclosed. 

2.  Integrated  Structural  Tailoring  and  Adaptive  Control  for  Advanced  Aircraft.  Wings,  L.  Li¬ 
brescu,  L.  Meirovitch  and  0.  Song* 

This  paper  presents  an  integrated  approach  combining  structural  tailoring  with  the  con¬ 
verse  piezoelectric  effect  for  the  purpose  of  actively  controlling  the  vibration  and  static 
aeroelastic  characteristics  of  advanced  aircraft  wings.  The  structural  model  incorporating 
a  number  of  nonclassical  features  consists  of  a  thin/thick- walled  closed  cross  section  can¬ 
tilevered  beam  whose  constituent  layers  exhibit  elastic  anisotropic  properties.  In  addition, 

*  The  authors  are  listed  in  alphabetical  order. 
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a  system  of  piezoelectric  actuators  bonded  to,  or  embedded  into  the  structure  generates  a 
localized  strain  field  in  response  to  an  injected  electric  current,  thus  producing  a  change  in 
the  static  and  dynamic  characteristics  of  the  structure.  Results  reveal  that  a  combination 
of  both  techniques  can  play  a  major  role  in  enhancing  the  vibrational  and  static  aeroelastic 
response  characteristics  of  aircraft  wings. 

The  paper  has  been  presented  at  the  AIAA/ ASME/ ASCE/ AHS/ ACS  34th  Structures, 
Structural  Dynamics,  and  Materials  Conference,  La  Jolla,  CA,  April  19-22,  1993,  and  is 
being  considered  for  publication  in  Journal  of  Aircraft.  A  copy  of  the  paper  is  enclosed. 

3.  Structural  Modeling  of  Low-Aspect  Ratio  Composite  Wings,  L.  Meirovitch  and  T.J.  Seitz 

This  paper  is  concerned  with  the  aeroelastic  tailoring  of  a  structural  model  consisting 

of  a  rigid  fuselage  and  a  low-aspect  ratio  wing  made  of  composite  materials.  The  wing  is 
modeled  as  a  trapezoidal  plate  with  root  and  tip  chords  parallel  to  the  flow  and  with  general 
sweep.  The  fuselage  is  capable  of  plunge  and  pitch  and  the  elastic  wing  model  includes  shear 
deformations  but  ignores  rotatory  inertia. 

The  paper  has  been  presented  at  the  AIAA/ASME/ASCE/ AHS/ACS  34th  Structures, 
Structural  Dynamics,  and  Materials  Conference,  La  Jolla,  CA,  April  19-22,  1993,  and  has 
been  tentatively  accepted  for  publication  in  Journal  of  Aircraft.  A  copy  of  the  paper  is 
enclosed. 

4.  Integrated  Structural  Design  and  Vibration  Suppression  Using  Independent  Modal  Space 

Control,  R.  A.  Canfield  and  L.  Meirovitch 

The  integrated  design  of  a  structure  and  its  control  system  is  treated  as  a  multiobjective 
optimization  problem.  Structural  mass  and  a  quadratic  performance  index  constitute  the 
vector  objective  function.  The  closed-loop  performance  index  is  taken  as  the  time  integral 
of  the  Hamiltonian.  Constraints  on  natural  frequencies,  closed-loop  damping,  and  actuator 
forces  are  also  considered.  Derivatives  of  the  objective  and  constraint  functions  with  respect 
to  structural  and  control  design  variables  are  derived  for  a  finite  element  beam  model  of 
the  structure  and  constant  feedback  gains  determined  by  independent  modal  space  control. 
Pareto  optimal  designs  generated  for  a  simple  beam  demonstrate  the  benefit  of  solving  the 
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integrated  structural  and  control  optimization  problem. 

The  .paper  has  been  presented  at  the  AIAA/ASME/ ASCE/AHS/ ACS  34^h  Strvctures, 
Structural  Dynamics,  and  Materials  Conference,  La  Jolla,  CA,  April  19-22,  1993,  and  has 
appeared  in  AIAA  Journal,  Vol.  32,  No.  10,  1994,  pp.  2053-2060.  A  copy  of  the  paper  is 
enclosed. 

5.  Vibration  and  Static  Aeroelastic  Instability  of  Nonuniform,  Thin-Walled  Beam  Compos¬ 
ite  Wings,  L.  Librescu,  L.  Meirovitch  and  0.  Song* 

The  equations  of  motion  for  a  nonuniform,  anisotropic  thin-walled  beam  are  derived 
and  applied  to  the  study  of  vibration  and  static  aeroelastic  instability  of  slender  tapered 
aircraft  wings  made  of  advanced  composite  materials.  Numerical  results  illustrate  the  effects 
of  anisotropy,  transverse  shear  flexibility,  primary  and  secondary  warping,  as  well  as  of  wing 
taper  ratio,  and  the  imphcations  of  these  effects  on  the  vibrational  and  divergence  instability 
characteristics  are  discussed. 

The  paper  has  been  presented  at  the  AIAA/ASME/ ASCE/AHS/ ACS  35th  Structures, 
Structural  Dynamics,  and  Materials  Conference,  Hilton  Head,  SC,  April  16-20,  1994,  and  is 
currently  being  revised  for  submission  for  publication  in  AIAA  Journal.  A  copy  of  the  paper 
is  enclosed. 


*  The  authors  are  listed  in  alphabetical  order. 
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Refined  Structural  Modeling  For  Enhanced  Vibrational 
and  Aeroelastic  Characteristics  of  Composite  Aircraft  Wings 

L.  Librescu,  L.  Meirovitch  and  0.  Song* 
Department  of  Engineering  Science  and  Mechanics,  Virginia  Polytechnic  Institute 

and  State  University,  Blacksburg,  VA  24061,  USA 

Abstract  An  analytical  study  focusing  on  the  vibrational  and  static  aeroelastic  response  of 
anisotropic  composite  aircraft  wings  modeled  as  thin-walled  beams  is  presented.  The 
theoretical  analysis  and  numerical  results  encompass  effects  such  as  transverse  shear 
flexibility  exhibited  by  the  advanced  composite  materials,  warping  restraint  charac¬ 
terizing  cantilevered  structures,  elastic  anisotropy  and  induced  structural  couplings. 
The  complex  effects  of  these  factors  are  highlighted  and  the  power  of  the  tailor¬ 
ing  technique  toward  enhancing  the  dynamic  and  static  structural  characteristics  is 
demonstrated. 

Keywords:  Composite  materials  -  Vibration  -  Static  aeroelasticity  -  Tailoring  - 
Transverse  shear  flexibility  -  Warping  restraint  -  Structural  couplings 

Resume 


*  The  authors  are  listed  in  alphabetical  order. 
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I  -  INTRODUCTION 


With  the  advent  of  high-performance  composite  materials,  thick-  and  thin-walled  beam 
structures  made  of  fiber  reinforced  laminated  composites  are  likely  to  play  an  increasing  role 
in  the  design  of  aircraft  wings.  In  addition,  the  ability  to  tailor  the  structural  characteristics 
argues  for  extensive  use  of  composite  materials  in  these  structures.  The  various  elastic 
and  structural  couplings  resulting  from  directional  properties  of  fiber-reinforced  composite 
materials  and  of  ply-stacking  sequence  can  be  exploited  successfully  to  enhance  the  response 
characteristics  of  flight  vehicles.  However,  extensive  use  of  these  potential  benefits  can  be 
achieved  only  if  the  effects  induced  by  these  couplings  are  well  understood. 

It  is  clear  that,  for  an  accurate  prediction  of  flight  vehicle  response  characteristics, 
comprehensive  structural  models  must  be  developed.  It  is  also  clear  that,  in  addition  to  this 
requirement,  more  powerful  analytical  tools  than  those  currently  available  are  needed  for 
accurate  prediction  of  structural  response  under  complex  static  and  dynamic  excitations. 

The  present  paper  is  devoted  to  the  development  of  the  structural  theories  and  ana¬ 
lytical  techniques  capable  of  treating  the  above  problem.  To  this  end,  a  refined  dynamical 
model  of  laminated  composite  thin-/thick- walled  beams  of  arbitrary  closed  cross  section  in¬ 
corporating  a  number  of  nonclassical  features  is  developed.  In  this  connection,  a  ply-angle 
scheme  generating  the  most  favorable  structural  coupling  for  the  problem  at  hand  is  imple¬ 
mented.  Then,  the  vibrational  and  static  aeroelastic  behavior  of  composite  swept-wings  is 
studied  and  the  influence  of  various  factors  on  the  response  characteristics  is  demonstrated. 

Two  methods  of  solution  for  the  problem  at  hand  were  explored  and  proved  to  be 
extremely  efficient.  An  exact  one,  based  on  the  Laplace  transform  technique  in  the  spatial 
domain,  and  an  approximate  one,  referred  to  as  the  extended  Galerkin  method,  yielding 
numerical  results  in  excellent  agreement  with  the  first.  The  analytical  developments  in  this 
paper  are  general  in  the  sense  that  they  are  valid  for  arbitrary  beam  cross  sections.  However, 
for  illustration  purposes,  a  biconvex  profile  typical  of  supersonic  wing  airplanes  is  adopted. 

It  should  be  remarked  that,  due  to  the  incorporation  of  transverse  shear,  the  ensuing 
beam  theory  is  perfectly  applicable  to  both  thin-  and  thick-walled  beams,  in  the  sense  of 
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^mai/^^0.1,  respectively,  where  hmax  denotes  the  maximum  thickness  of  the  wall  and  h  a 
typical- cross-sectional  dimension.  Although  the  theory  is  applicable  to  thick-walled  beams 
as  well,  the  standard  terminology  of  thin- walled  beams  is  used  throughout.  It  should  be 
stressed  here  that  the  necessity  of  incorporating  transverse  sheer  effects  arises  not  only  from 
the  fact  that  composite  beams  tend  to  be  thicker  than  the  standard  metallic  counterpart, 
but  also  from  the  fact  that  advanced  fiber  composite  materials  exhibit  high  flexibilities  in 
transverse  shear. 

It  should  be  emphasized  here  that  the  model  of  composite  aircraft  wings  almost  uni¬ 
versally  used  to  date  has  been  the  simple  solid  beam  [1].  Hence,  the  composite  thin- walled 
beam  model  developed  and  used  in  this  paper  represents  a  significant  advance  in  the  state  of 
the  art  as  it  permits  not  only  the  incorporation  of  effects  not  accounted  for  previously,  but 
also  enables  one  to  treat  a  number  of  important  problems  of  wing  structures  using  a  more 
refined  structural  model. 

II  -  BASIC  ASSUMPTIONS 

This  study  uses  the  structural  concept  of  single-cell,  thin-walled  beams  of  arbitrary 
cross-section.  Pertinent  quantities  are  referred  to  a  global  coordinate  system  xyz,  where  x 
and  y  denote  the  cross-sectional  coordinates  and  z  is  the  spanwise  coordinate  (Fig.  la).  The 
theory  is  based  on  the  following  assumptions;  i)  beam  cross-sections  do  not  deform  in  their 
own  planes,  ii)  transverse  shear  effects  are  significant,  iii)  the  twist  varies  along  the  span,  i.e., 
the  rate  of  twist  d<i>ldz  is  no  longer  assumed  to  be  constant  (as  in  the  Saint- Venant  torsional 
model)  but  a  function  of  the  spanwise  coordinate,  where  (j)'  —  d(l>ldz  constitutes  a  measure 
of  the  torsion-related  nonuniform  warping,  iv)  primary  and  secondary  warping  effects  are 
sufficiently  important  to  be  included.  (The  first  one  is  related  to  the  warping  displacement 
of  points  on  the  midline  cross  section  and  the  second  one  is  related  to  points  off  the  midline 
contour.)  and  v)  in  the  absense  of  an  internal  pressure  field,  the  hoop  stress  resultant  is 
negligibly  small  compared  to  the  remaining  stresses. 

III  -  KINEMATICAL  EQUATIONS 

Consistent  with  assumption  iv),  the  primary  warping  function  Fu,{s)  is  taken  as  [2-4] 
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where 


,  §rn{s)ds  2Ac 

denotes  the  torsional  function,  Ac  the  cross-sectional  area  of  the  beam  bounded  by  the 
midline  contour,  s  the  arc-length  measured  along  the  circumferential  coordinate,  whose  ori¬ 
gin  is  arbitrarily  but  conveniently  chosen,  s  is  a  dummy  coordinate  associated  with  the 
s-coordinate,  f(-)ds  the  integral  along  the  closed  midline  contour,  r„(s)  =  x(s)^  +  y{s)m  a 
geometric  quantity  (see  Fig.  lb),  where  £  =  cos(n,  x),  m  =  cos(n,  y)  denote  direction  cosines, 
and  /?  =  /  ds  is  the  circumferential  contour  length. 

In  accordance  with  assumptions  i)  -  iv),  and  in  order  to  reduce  the  three-dimensional 
elasticity  theory  of  thin-walled  beams  to  an  equivalent  one- dimensional  one,  the  components 
of  the  displacement  vector  are  expressed  as  [4] 

u{x,  y,  2,  t)  =Uo{z,  t)  -  y<l){z,  t),  v(x,  y,  z,  t)  =  Vo(z,  t)  +  x(^(z,  t)  (4a,  b) 

w{x,  y,  2,  t)  =Wo{z,  t)  -h  6x{z,  t)  [y(5)  -f-  nm] 

-f  Oy{z,  t)  [x(s)  -I-  n£]  —  <j){z,  t)  [Fu;(s)  -f  na(s)]  (4c) 

In  the  above  equation,  and  na{s)  play  the  role  of  primary  and  secondary  warping  func¬ 

tions  (see  also  [5]),  respectively,  where  n  denotes  the  coordinate  in  the  thickness  direction, 
and 


=  lyz{z,i)-Vc{z,t),  ey{z,t)  =  ^xz{z,t)-u^{z,t),  a{s)  = -y{s)£  +  x{s)m  (5a-c) 

in  which  di(z,t)  and  Oy{z,t)  denote  the  rotations  about  axes  x  and  y,  respectively,  and 
jyz  and  7x2  are  the  transverse  shear  strains  in  the  planes  yz  and  xz,  respectively;  we  note 
that  primes  denote  derivatives  with  respect  to  z.  When  transverse  shear  effects  are  ignored. 
Ox  — >  —Vq  and  Oy  —*■  —u'q.  In  agreement  with  Eqs.  (4),  it  can  be  readily  verified  that  the 
assumption  of  cross-sectional  nondeformability,  implying  Cxx  =  0,  tyy  =  0,  and  txy  =  0,  or 


equivalently  e„„  =  633  =  ejn  =  0,  and  the  continuity  requirement  of  w  along  the  midline 
contour,  i.e.,  §{dwfds)ds  =  0,  are  fulfilled. 

The  quantities  UQ{z,t),  vo{z,t)  and  wo{z,t),  denoting  the  rigid-body  translations  along 
the  x,y  and  2:  axis,  respectively,  and  0x{z,t),6y{z,t)  and  <j){z,t)^  denoting  the  rigid-body 
rotations  about  the  x-  and  y-axes  and  the  twist  about  the  2-axis,  respectively,  constitute  the 
variables  of  the  problem. 


IV  -  CONSTITUTIVE  EQUATIONS 


Consider  the  case  of  composite  thin-walled  beams  consisting  of  a  finite  number  N  of  ho¬ 
mogeneous  layers.  It  is  assumed  that  the  material  of  each  constituent  layer  is  linearly  elastic 
and  anisotropic  and  that  the  bonding  between  the  layers  is  perfect.  The  three-dimensional 
constitutive  equations  for  a  generally  orthotropic  elastic  material  can  be  expressed  as 
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where  Qij  denote  the  transformed  elastic  coefficients  associated  with  the  A:th  layer  in  the 
global  coordinate  system  of  the  structure  and  7pr  =  2tpr,  p  ^  r  and  e,y  denote  the  com¬ 
ponents  of  the  strain  tensor.  The  three-dimensional  dependence  in  Eqs.  (6)  can  be  re¬ 
duced  to  an  equivalent  one-dimensional  dependence  in  two  steps.  The  first  step,  yielding 
the  two-dimensional  local  constitutive  equations,  consists  of  the  integration  of  the  original 
three-dimensional  form  through  the  laminate  thickness,  while  the  second  step,  resulting  in 
the  one-dimensional  form,  consists  of  the  integration  of  the  previous  form  of  constitutive 
equations  along  the  midline  contour  of  the  beam  cross  section. 

In  light  of  assumptions  i)  and  v),  the  local  constitutive  equations,  expressed  in  terms 
of  the  strain  measures,  become  the  equations  for  the  membrane  stress  resultants 


Nzz  =  Ki\e^xz  +  Kuisz  + 

N3Z  =  Kiielx  +  K22isz  +  K2z<f>'  +  K24e:x 
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(7a) 

(76) 


transverse  shear  stress  resultant 


^zn  —  ^A47zn 

(8) 

and  stress  couples 

Lzz  =  -f  K^2laz  +  Ka3<I>  +  A44e”^ 

(9a) 

Lsz  =  -t-  Ks27sz  +  +  KsaF^z 

(96) 

where  Kij  denote  the  modified  local  stiffness  coefficients,  (see  Appendix).  Consistent  with 
Eqs.  (4)  and  (5),  the  strain  components  entering  into  the  above  constitutive  equations  are 
Czz  =  +  ne”^,  ')sz  =  Isz  +  Isz  and  7^„,  where 

^°zz  =^o  +  +  p(s)d'^  -  F^(s)<^",  +  0'^m  -  a(s)<f>"  (10a,  b) 

isz  =  -  (n'o  +  Oy)m  +  (vq  +  7sz  =  O^zn  =  (n'o  +  &y)£  +  (uq  +  0i)m  (lOc-e) 

In  the  above  equations,  and  e”^  denote  the  axial  strain  components  associated  with 
the  primary  and  secondary  warping,  respectively,  and  stand  for  the  tangential  shear 
strains  of  the  beam  midsurface  induced  by  transverse  shear  and  by  twist,  respectively,  and  ')zn 
denotes  the  transverse  shear  strain  component.  The  stress  resultants,  the  stress  couples  and 
the  strain  measures  appearing  in  Eqs.  (7)  -  (9)  exhibit  a  two-dimensional  spatial  dependence, 
the  dependence  being  on  s-  and  2-coordinates.  In  the  dynamical  problem,  they  also  depend 
on  time. 


V  -  THE  BOUNDARY- VALUE  PROBLEM 

The  boundary- value  problem,  consisting  of  the  differential  equations  and  the  boundary 
conditions,  can  be  derived  conveniently  by  means  of  the  extended  Hamilton’s  principle,  which 
can  be  stated  as  follows  [6]; 

/  (6T  —  SV  -h  8W)  dt  =  0,  Suq  =  Svq  =  Swq  =  60x  =  dOy  =  8(f)  =  0  aX  t  =  ti,t2  (11) 

Jti 

where  T  and  V  denote  the  kinetic  energy  and  strain  energy,  respectively,  while  8W  is  the 

virtual  work  due  to  nonconservative  forces.  The  kinetic  energy  has  the  form 
N 


1  ^  r 

-  0  I  {uo  -  yf)'^  +  [vo  +  +  \^(wo  +  yOx  +  xOy  - 

-|-n  (iOy  -|-  mdx  —  I  dndsdz 


(12) 
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whereas  the  strain  energy  can  be  shown  to  have  the  expression 

+y6'j.  -  Fu(j>"  +  n  {l6'y  +  mO'j.  -  -  ^Uq  +  Oy)  m 

+  (y'o  +  ^i)  ^  [(“o  +  ^y)  ^  +  (^0  +  I  dndsdz 


(13) 


where  L  denotes  the  wing  semispan.  Moreover,  the  virtual  work  of  the  nonconservative  forces 
can  be  written  as 


SWz= 


+  mzS(j>)  dz 


where  py  and  denote  the  hft  force  per  unit  length  and  aerodynamic  twist  moment  (positive 
nose  up)  about  the  elastic  axis. 

Carrying  out  integrations  with  respect  to  n,  s  and  t,  we  can  write 


[  STdt  =  -f  I  /  {IiSuo  +  I26V0  +  hSwo 
Jto  Jto  [v/o 

+  (h  - 19)  H  +  hSOy  +  hse^)  dz  -  i9<f> 


dt 


(15) 


where  /,  denote  the  inertia  terms,  given  in  the  Appendix.  Integration  with  respect  to  n  and 
s  in  Eq.  (13)  and  consideration  of  Eqs.  (1),  (5),  (7)  and  (10),  yields 

8V  =  -  [r.SwQ  +  [M'y  -  Q^)  SOy  +  (m'  -  Qy)  69,  +  +  M')  6<l> 

J  0 

dz  +  +  MySOy  +  MxdO,  —  B^^6<j) 

+  6(f)  +  Qxduf)  +  (16) 


where 

Tz{z,t)  =  Nzzds,  Qx{z,t)  =  j>  {-Nszm  +  Nzn^)ds 
Qy{z,t)=  <f  {Nsz^  +  Nznm)ds,  Mx{z,t)=  {yNzz  +  Lzzm)ds 

\  /  (17a  -  9) 

My{z,  t)  =  {xNzz  +  Lzzt} ,  Mz{z,  t)  =  2  j)  Nazifjds 
Bc^{z,t)  =  j  [F^{s)Nzz  +  a{s)Lzz]ds 

In  Eqs.  (17), Tz,Qx  and  Qy  denote  the  axial  force  and  shear  forces  in  the  x-  and  y-  directions, 
Mx,  My  and  Mz  denote  the  moments  about  the  x-,  y-  and  2:-axes,  respectively,  and  B^ 
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denotes  the  bimoment.  Introducing  Eqs.  (14)-(16)  into  Eq.  (11)  and  following  the  usual 
steps  [I]  ,  we  obtain  the  boundary- value  problem  for  the  most  general  case  of  anisotropy.  The 
boundary- value  problem  consists  of  six  differential  equations  of  motion  for  the  displacements 
uo,vo,WoiOx.,Oy  and  4>,  together  with  the  corresponding  boundary  conditions.  Such  a  set 
exhibits  complete  coupling  between  the  various  modes,  i.e.,  primary  and  secondary  warping, 
vertical  and  lateral  bending,  twist  and  transverse  shear.  However,  the  principal  goal  of 
structural  tailoring  lies  in  the  appropriate  selection  of  fiber  orientation  so  as  to  produce 
desired  elastic  couplings  between  certain  modes.  For  the  problem  at  hand,  the  induced 
elastic  couplings  must  play  a  decisive  role  in  the  enhancement  of  the  free  vibration  and 
aeroelastic  response  characteristics  of  wing  structures.  At  the  same  time,  the  selected  ply- 
angle  orientation  should  not  generate  undesired  couplings,  producing  unwanted  effects  on 
the  response  characteristics.  In  this  sense,  the  bending-twist  coupling  is  most  important  in 
the  design  of  aircraft  wings.  Its  beneficial  effects,  demonstrated  for  the  solid  beam  model 
in  [1],  will  be  considered  here  in  the  context  of  a  thin- walled  beam  model.  Additional 
beneficial  effects  of  this  cross- coupling  have  been  highlighted  recently  in  [8].  The  above 
criteria  for  selecting  fiber  orientation,  together  with  ease  of  implementation  in  design  and 
manufacturing,  result  in  the  ply-angle  distribution 

^y)  =  -H-y)-  (18) 

This  ply  configuration  is  shown  in  Fig.  2.  According  to  terminology  used  in  [9],  struc¬ 
tures  displaying  this  ply-angle  distribution  are  said  to  exhibit  circumferentially  asymmetric 
stiffness  configuration. 

In  the  case  of  the  ply- angle  distribution  given  by  Eq.  (18),  the  extended  Hamilton’s 
principle,  Eq.  (11),  yields  two  independent  boundary- value  problems.  The  first  boundary- 
value  problem  is  of  eighth  order  and  involves  the  coupling  of  the  twist  <f,  the  vertical  bending 
Vo  and  the  fiapwise  transverse  shear  6x.  On  the  other  hand,  the  second  boundary- value  is 
of  sixth  order  and  involves  the  coupling  of  the  extension  luo,  the  lateral  bending  uq  and 
the  chordwise  transverse  shear  6y.  The  first  boundary- value  problem  is  governed  by  the 
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differential  equations  of  motion 


-aee(f>""  +  ajsO"  -  a65(^>o'  +  O'x)  +  +  ruz  =  (64  +  65)^  -  (610  +  hs)’^"  (19a) 

assiv'o  +  O'j.)  +  056^  + Py  =  bivo  (196) 

033^1  +  Ci37<f>"  -  aesiv"  +  O'x)  -  as6<t>"  =  (64  +  h^)6x  (19c) 

to  be  satisfied  over  ^  <  z  <  L.  For  cantilevered  thin-walled  beams,  the  solution  of  Eqs.  (19) 
must  satisfy  the  boundary  conditions 

(^  =  0,  no  =  0,  ^1  =  0,  ^  =  0  (20a-d) 

at  z  =  0  and 

—  +  073^1  —  Q65(t^0  +  O'x)  +  =  —(610  -f  6ig)^' 

(21a-d) 

055(^^o  +  Ox)  -b  a^(,<f>"  —  0,  033^^  -b  =  0,  aee(l>"  +  055(^0  +  Ox)  =  0 
at  z  =  i.  Note  that  singly  and  doubly  underlined  terms  in  Eqs.  (19)-(21)  are  associated 
with  the  warping  restraint  and  warping  inertia  effects,  respectively. 

The  couplings  exhibited  by  the  second  boundary- value  problem  are  of  no  interest  here, 
so  that  the  problem  will  not  be  pursued  any  further. 

The  boundary-value  problem  associated  with  the  bending-twist-transverse  shear  mo¬ 
tions,  Eqs.  (19)-(21),  can  be  used  to  enhance  the  vibrational  and  aeroelastic  response  char¬ 
acteristics  of  wing  structures.  It  should  be  mentioned  here  that  the  stiffness  terms  037  =  073 
and  056  =  065  appearing  in  Eqs.  (19)  and  (21)  are  responsible  for  the  coupling  between 
bending  and  twist,  with  the  effect  of  055  due  to  warping  being  much  weaker  than  the  effect 
of  037  (see  Fig.  11).  It  should  be  pointed  out  that  deletion  of  the  warping  restraint  results  in 
a  reduction  in  the  order  of  the  boundary-value  problem  from  eight  to  six,  so  that  only  three 
boundary  conditions  must  be  satisfied  at  each  end.  On  the  other  hand,  if  transverse  shear 
effects  are  ignored,  then  the  order  of  the  boundary- value  problem  is  preserved.  The  stiffness 
coefficients  a,j  and  the  inertia  coefficients  6,-  appearing  in  Eqs.  (19)  and  (21)  are  displayed 
in  the  Appendix. 
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VI  -  STRUCTURAL  TAILORING  FOR  IMPROVED  VIBRATION  AND 
STATIC  AEROELASTIC  RESPONSE 

Static  aeroelastic  behavior,  which  includes  both  divergence  instability  and  aeroelastic 
lift  distribution  is  an  important  consideration  in  the  design  of  modern  aircraft.  The  analysis 
performed  here  addresses  the  problem  of  designing  the  wing  so  as  to  take  advantage  of 
structural  couplings  from  a  static  aeroelastic  viewpoint.  This  is  done  by  using  the  unique 
directional  properties  of  advanced  composite  materials.  The  same  importance  should  be 
afforded  to  the  vibrational  characteristics,  which  are  basic  to  the  dynamic  response,  flutter 
instability  and  aeroservoelasticity  studies  of  flight  vehicles. 

In  the  case  of  free  vibration,  the  terms  associated  with  the  external  loadings  are  omit¬ 
ted.  Moreover,  for  static  aeroelastic  problems,  the  inertia  terms  must  be  discarded  from 
Eqs.  (19)  and  (21),  and  the  only  loading  terms  to  be  retained  are  the  ones  associated 
with  the  aerodynamic  lift  py  and  the  torsional  aerodynamic  moment  m^.  Using  strip-theory 
aerodynamics,  we  can  write  [10] 

Py{z)  =qncao{(j)o  +  (j)-  v'g  tan  A)  —  NWf2L  (22a) 

=qnCaoe{<f)o  +  <P  -  v'g  tan  A)  +  qnC^CMAC  -  NWdl2L  (226) 

Here  q^  =  \pU^  denotes  the  dynamic  pressure  normal  to  the  leading  edge  of  the  swept  wing, 
c  the  chord  of  the  wing,  Ug  the  corrected  lift  curve  slope  coefficient,  A  the  angle  of  sweep 
(considered  positive  for  swept-back),  e  the  offset  between  the  aerodynamic  and  reference 
axis,  (j>g  the  rigid  angle  of  attack  (measured  in  planes  normal  to  the  leading  edge).  Cm  AC  the 
wing  section  pitching  moment  coefficient  (whose  influence,  as  usual,  is  disregarded),  VU/2Z 
the  wing  weight  per  unit  length  and  N  the  load  factor  normal  to  the  wing  surface,  whose 
expression  is 

N  =  f  {(jfg  +  (l>  —  v'gt&n  /i)dz  (23) 

W  Jo 

The  static  aeroelastic  response  is  analyzed  here  both  in  the  subcritical  range,  i.e., 
for  the  range  of  velocities  <  (^n)/?,  where  {qn)D  denotes  divergence  dynamic  pressure, 
and  in  the  critical  case  as  well.  As  a  general  remark,  Eqs.  (22)  reveal  that  for  A  <  0, 
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i.e.,  for  swept-forward  wings,  the  aeroelastic  bending-twist  coupling  results  in  an  increase 
in  py(0)  and  mz{z)^  which  in  turn  reduces  dramatically  the  divergence  speed,  whereas  for 
A  >  0,  i.e.,  for  swept-back  wings,  the  opposite  trend  takes  place.  These  two  phenomena  are 
referred  to  as  wash-in  and  wash-out  [11],  respectively.  Whereas  the  goal  of  the  subcritical 
aeroelastic  analysis  consists  of  the  determination  of  the  distribution  of  the  effective  angle 
of  attack  <^e//  of  fh®  lift  force,  as  affected  by  the  elastic  deformations,  the  study  of 
the  critical  case  involves  the  determination  of  divergence  instability  conditions.  Clearly,  the 
main  target  of  tailoring  applied  to  swept-forward  wing  is  to  yield  a  decrease  of  the  effective 
angle  of  attack,  and  implicitly  of  the  aeroelastic  lift  and,  as  a  result,  an  increase  of  the 
critical  divergence  speed.  Whereas  the  study  of  the  subcritical  static  aeroelastic  response 
requires  the  solution  of  an  integral-differential  system  of  equations,  obtained  by  inserting 
Eqs.  (22)  and  (23)  into  Eqs.  (19),  the  determination  of  the  divergence  speed  leads  to  the 
solution  of  an  eigenvalue  problem,  where  the  divergence  speed  plays  the  role  of  eigenvalue. 
Structural  tailoring  applied  to  the  vibration  of  wing  structures  must  result  in  an  increase 
in  the  eigenfrequencies  without  weight  penalties.  To  determine  the  natural  frequencies,  one 
must  solve  an  eigenvalue  problem. 

In  spite  of  the  mathematical  complexities  involved,  the  two  previously  mentioned  solu¬ 
tion  techniques  used  here  proved  to  be  extremely  powerful,  as  demonstrated  by  the  following 
numerical  illustrations.  Details  of  the  techniques  can  be  found  in  [12-14]. 

VII  -  SPECIAL  CASES  INVOLVING  DIVERGENCE  INSTABILITY 

In  the  most  general  case,  closed-form  solutions  for  the  divergence  speed  are  not  feasible. 
However,  in  a  number  of  special  cases  closed-form  solutions  can  be  obtained.  These  cases 
are  concerned  with  i)  pure  bending  divergence  of  swept  wings  infinitely  rigid  in  transverse 
shear  and  ii)  pure  torsion  divergence. 

In  the  first  case,  eliminating  from  Eqs.  (196)  and  (19c),  eissuming  very  large 

torsional  stiffness,  letting  Ox  — >•  —Vq,  ignoring  the  inertia  terms  and  implementing  a  Rayleigh- 
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quotient  procedure,  the  expression  of  the  divergence  speed  can  be  shown  to  have  the  form 

2033  /  (vq) 

(qn)D  =  - 


^dz 


oqc tan  A  /  {vl)'dz 

Jo 

where  vq  must  satisfy  the  boundary  conditions 

uq  =  0,  Uq  =  0  at  2:  =  0 
Vq  =  0,  Vq'  =  0  a.t  z  =  L 


(24) 


(25a) 

(256) 


Equation  (24)  reveals  that  only  swept-forward  wings,  A  — >  —A,  can  exhibit  divergence 

instability  in  pure  bending.  This  result  represents  an  extension  to  wings  modeled  as  thin- 

walled  beams  of  results  obtained  in  [16].  Also  from  Eq.  (24),  one  can  conclude  that  the 
bending  stiffness  term  033  must  be  maximized  to  increase  {qn)D  ^.s  much  as  possible.  This 
can  be  achieved  by  designing  the  wing  so  that  Ku  — >  0. 

In  the  pure  torsion  case,  assuming  infinite  bending  stiffness  and  implementing  a 
Rayleigh-quotient  procedure  in  conjunction  with  the  static  counterpart  of  Eq.  (19a),  it  can 

be  shown  that  ^ 

I  \a66{<l>"f  +  dz 

{qn)D  =qD=  L  (26) 

ooce  /  (jP'dz 

Jo 

in  which  the  underdashed  term  is  connected  with  the  warping  inhibition,  and  note  that  (f) 
must  satisfy  the  boundary  conditions 


</>  =  0,  <f)'  =  0  ai  z  =  0 

<j)"  =  a^^(t)  —  aQQ^'"  —  0  &t  z  =  L 


(27a) 

(276) 


As  in  the  case  of  wings  modeled  a.s  solid  beams  [16],  Eq.  (26)  reveals  that  pure  torsion 
divergence  can  occur  for  straight  wings  only. 

VIII  -  APPROXIMATE  EXPRESSION  FOR  THE  DIVERGENCE  OF  SWEPT- 
FORWARD  WINGS 

An  approximate  expression  for  the  coupled  divergence  of  swept-forward  wings  can  be 
derived  using  the  conditions  corresponding  to  decoupled  divergence  in  bending  and  torsion. 
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Eqs.  (24)  and  (26).  The  expression  is  based  on  a  linear  algebraic  relation  between  the  two 
decoupled  expressions  of  the  divergence  instability  obtained  previously.  The  linear  relation 
yields  the  expression  for  the  divergence  in  the  form 


Md  = 


1  P 
S 


(28) 


where 


where  rj  =  z/L.  As  in  Eqs.  (24)  and  (26),  the  warping  restraint  effect  is  included  but 
the  transverse  shear  flexibility  is  ignored.  Equation  (28)  represents  the  extension  to  wings 
modeled  as  thin- walled  beams  of  the  divergence  expression  obtained  in  [16]  for  solid  beams. 


IX  -  NUMERICAL  ILLUSTRATIONS 

Wing  structures  modeled  as  cantilevered  thin- walled  beams  of  biconvex  cross  sections 
made  of  the  graphite-epoxy  material  are  considered.  The  material  properties  used  are 


El  =30  X  10®  psi,  Et  =  0.75  x  10®  psi 
Git  =0.37  x  10®  psi,  Grr  =  0.45  x  10®  psi 
y^LT  —yiT  =  0.25,  p  =  14.3  X  10”®  lb  sec^  /in^ 

where  L  and  T  denote  directions  parallel  and  transverse  to  the  fibers,  respectively. 

The  geometrical  characteristics  of  the  wing  are  displayed  in  Fig.  la.  Figure  3a  shows 
the  first  three  in  vacuo  eigenfrequencies  associated  with  Problem  A  as  functions  of  the  ply 
angle  6  with  the  bending-twist  coupling  stiffness  first  included,  037  ^  0,  and  then  discarded, 
037  =  0.  In  both  cases,  plots  of  the  frequencies  versus  6  are  symmetric  about  6  =  90° 
and  experience  as  many  peaks  as  the  eigenfrequency  number.  When  037  =  0,  the  second 
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eigenfrequency  has  a  local  maximum  in  the  vicinity  of  9  =  75°  and  6  =  105°,  and  at  these  ply 
angles  it  comes  in  close  proximity  to  the  third  eigenfrequency.  On  the  other  hand,  when  the 
effect  of  bending-twist  coupling  is  included,  implying  037  ^  0,  Fig.  3a  reveals  that  frequency 
near  merging  is  precluded.  This  phenomenon  was  also  observed  in  Ref.  15,  in  which  a  solid 
beam  model  was  used.  The  frequency  of  the  fundamental  mode  has  a  maximum  at  0  =  90°. 
At  0  =  0°  (and  180°),  where  the  bending  and  twist  become  decoupled,  the  fundamental 
mode  can  be  identified  as  a  pure  bending  mode,  denoted  by  B.  The  frequency  associated 
with  the  second  mode  first  increases  for  0°  <  0  <  80°  and  then  decreases  for  80°  <  6  <  90°, 
the  trend  being  symmetric  about  6  =  90°.  Another  notable  trend  is  that,  in  the  absence  of 
cross- coupling  rigidity,  the  second  eigenfrequency  is  overestimated  compared  to  the  real  case, 
in  which  the  cross-coupling  rigidity  is  included.  At  0  =  0°  (or  9  =  180°)  and  9  =  90°,  where 
decoupling  occurs,  the  second  mode  can  be  identified  as  the  second  pure  bending  mode  and 
the  first  pure  torsional  mode,  respectively,  where  the  latter  is  denoted  by  T.  As  far  as  the 
third  eigenfrequency  is  concerned,  a  more  complex  variation  with  the  ply  angle  is  observed. 
Indeed,  the  variation  about  9  =  90°  is  drastically  attenuated  when  the  cross  coupling  is 
ignored.  At  0  =  0°  (or  9  =  180°)  and  9  —  90°,  this  mode  can  be  identified  as  the  first  pure 
torsional  mode  and  second  pure  bending  mode,  respectively. 

Figure  3b  displays  the  first  three  eigenfrequencies  for  Problem  A  as  functions  of  the 
ply  angle  for  the  cases  in  which  the  transverse  shear  effects  are  incorporated  and  ignored. 
Clearly,  for  ply  angles  such  that  the  bending  is  dominant,  omission  of  the  transverse  shear 
causes  an  overestimation  of  frequencies. 

Figure  4a  portrays  plots  of  the  first  three  eigenfrequencies  for  Problem  B  versus  the 
ply  angle  9.  In  this  case,  the  lateral  bending-extension  cross-coupling  parameter  034  plays  a 
similar  role  to  the  bending-twist  stiffness  037.  However,  these  three  eigenfrequencies  are  well 
separated  within  the  entire  range  of  ply  angles  9  and,  in  addition,  removal  of  034  does  not 
cause  the  second  and  third  eigenfrequencies  to  approach  each  other,  regardless  of  the  value  of 
9.  The  plots  are  symmetric  about  9  =  90°  and,  as  expected,  the  lateral  eigenfrequencies  are 
much  higher  than  the  transverse  counterpart  frequencies.  At  0  =  90°  the  first  three  modes 
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are  identified  as  the  first  three  pure  lateral  bending  modes,  denoted  by  B,  and  at  0  =  0°  and 
9  =  180°  as  the  first  two  pure  lateral  bending  modes  and  the  first  pure  axial  mode,  where 
the  latter  is  denoted  by  A. 

Figure  4b  portrays  the  first  three  eigenfrequencies  in  lateral  vibration  versus  the  ply 
angle  9  for  the  cases  in  which  the  transverse  shear  is  included  and  ignored.  A  comparison  of 
Figs.  3b  and  4b  reveals  that  the  overestimation  of  eigenfrequencies  for  the  case  in  which  the 
transverse  shear  effects  are  ignored  is  more  pronounced  in  chordwise  vibration  than  in  the 
flapping  vibration  counterpart. 

Figure  5  presents  the  bending  and  twist  in  the  three  lowest  modes  as  functions  of  the 
normalized  position  t)  =  zf  L  for  9  =  45°.  The  modes  are  normalized  so  that  the  value  at  the 
tip  is  equal  to  unity.  For  other  ply  angles  and  for  the  case  of  lateral  bending,  similar  plots 
were  obtained,  but  they  are  not  displayed  here.  As  a  general  comment,  the  position  of  the 
nodal  points  changes  in  general  with  9.  This  trend,  coupled  with  large  variations  of  these 
modes  with  9,  is  likely  to  have  a  significant  effect  on  the  flutter  behavior. 

The  effect  of  the  ply  angle  on  the  divergence  speed  is  illustrated  in  Fig.  6  for  the  swept- 
back,  swept-forward  and  straight  wing  aircraft.  The  figure  shows  plots  of  the  normalized  di¬ 
vergence  speed  {qn)j}  /  (9n)^  versus  the  ply  angle  9,  where  (g„)^  =  (^n)^)  =  35.566psi. 

It  reveals  that  the  range  of  ply  angles  for  which  divergence  instability  is  avoided  decreases 
with  increasing  forward  sweep  and  increaises  with  increasing  sweep  back  angles.  Figure  7 
displays  the  normalized  divergence  speed  versus  the  ply  angle  for  various  values  of  the  wall 
thickness  and  for  A  =  0  and  A  =  —60°.  It  reveals  that,  when  an  increase  in  the  wall  thickness 
is  an  option  for  increasing  the  divergence  speed,  then  it  must  be  considered  in  conjunction 
with  the  proper  ply  angle,  thus  yielding  a  maximum  increase  of  the  divergence  speed. 

In  Fig.  8a  plots  of  the  divergence  speed  versus  the  ply  angle  for  two  swept  angles, 
A  =  0°  and  A  =  —60°,  and  for  the  cases  037  =  0  and  037  ^  0  are  displayed.  As  the  figure 
reveals,  for  037  =  0  the  divergence  speed  exhibits  symmetry  about  9  =  90°,  whereas  for 
037  ^  0  it  is  not  symmetric.  Moreover,  results  for  straight  wings  reveal  that  a  larger  range 
of  ply  angles  corresponding  to  infinite  divergence  speed  is  exhibited  than  for  the  forward- 
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swept  wing  counterpart.  Although  the  bending-twist  cross  coupling  is  not  the  only  factor 
influencing  the  aeroelastic  static  behavior,  it  is  clear  that  its  effect  is  dominant.  By  increasing 
its  negative  value  as  much  as  possible,  the  bending  of  the  swept  wings  and  the  resulting  twist 
tend  to  produce  a  wash-out  effect. 

Figure  8b  displays  a  plot  of  the  ratio  (?„) ^  /  (qn)])  WR  ^  function  of  the  ply  angle 

^  for  A  =  —30°  and  A  =  0,  where  {qn)D  pw  divergence  speed  for  a  free  warping  model 
and  {qn)DWR  divergence  speed  for  the  case  in  which  the  warping  restraint  effect  is 

included.  It  essentially  demonstrates  the  effect  of  warping  restraint  on  the  divergence  speed. 
It  can  be  concluded  from  this  figure  that  the  warping  inhibition  does  not  always  result  in  an 
increase  in  the  divergence  speed.  In  other  words,  the  free  warping  model  does  not  yield  the 
most  conservative  results  from  the  divergence  instability  point  of  view.  Hence,  for  a  rational 
design,  the  warping  restraint  effect  must  always  be  taken  into  consideration.  A  similar  result 
was  obtained  for  a  solid  beam  wing  model  in  [16-18]. 

In  view  the  importance  of  this  cross-coupling  rigidities,  in  Fig.  9  plots  of  037,  033,  077, 
056  and  014  versus  0  are  portrayed.  Based  on  these  plots,  a  value  of  6  resulting  in  a  more 
rational  design  can  be  achieved.  Within  the  theory  of  solid  beams,  the  importance  of  the 
bending-twist  coupling  was  underscored  in  [1,11-13,  15-17] 

Figures  10a  and  10b  display  plots  of  the  normalized  effective  angle  of  attack,  4>eff/<f>0, 
versus  the  normalized  position  t]  =  zjL  for  a,  swept  and  a  straight  wing,  and  with  the  ply 
angle  9  and  sweep  angle  A  acting  as  parameters,  where  =  (t>Q  +  <f>—VQ  tan  A  is  the  effective 
angle  of  attack  and  (f>Q  =  5°  is  the  rigid  angle  of  attack.  The  plots  provide  a  measure  of  the 
subcritical  static  aeroelastic  response.  The  effective  angle  of  attack  constitutes  a  measure 
of  the  induced  aeroelastic  loads.  For  ply  angles  9  <  90°,  the  aeroelastic  loads  are  amplified, 
whereas  for  9  >  90°,  for  both  forward  swept  and  swept  back  wings,  they  are  attenuated. 
In  all  cases,  the  flight  speed  corresponds  to  a  dynamic  pressure  (9n)flight  =  3.446  psi.  The 
basic  conclusion  is  that  tailoring  can  play  a  significant  role  not  only  in  counteracting  the 
detrimental  wash-in  effect,  but  also  in  diminishing  the  effect  of  elastic  twist.  This  change  of 
the  traditional  subcritical  static  aeroelastic  response  of  forward  swept  wing  is  basically  due 
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to  the  bending-twist  coupling  stiffness  037. 

The  results  obtained  here  confirm  that  the  coupling  stiffness  037  plays  a  key  role  in  con¬ 
trolling  the  wing  behavior  as  far  as  the  subcritical  static  aeroelastic  response  and  divergence 
instability  are  concerned.  By  increasing  the  negative  value  of  037  as  much  as  possible  for 
A  <  0,  the  wash-in  effect  turns  into  a  wash-out  effect.  The  effect  of  transverse  shear  flexibil¬ 
ity  on  the  subcritical  aeroelastic  response  of  swept-forward  wings  is  in  general  detrimental, 
in  the  sense  that  it  exacerbates  the  wash-in  effect.  This  conclusion,  based  on  results  obtained 
in  this  study  but  not  displayed  here,  is  reinforced  by  results  obtained  in  [13].  However,  in 
contrast  to  other  cases  of  ply  angles,  results  for  55°  <  6  <  90°  and  for  A  >  0  reveal  that 
the  coupling  rigidities  play  a  significant  role  also  in  this  respect  and  can  render  the  effect  of 
transverse  shear  flexibility  either  immaterial  or  slightly  beneficial. 

Finally,  Fig.  11  displays  the  subcritical  aeroelastic  responses  of  Vq  =  vq/L,  Ox  and 
(j)  as  functions  of  tj  for  a  straight  wing  with  ply  angles  6  =  45°  and  6  =  135°.  As  before, 
(9n)fiight  =  3.448  psi  was  prescribed.  These  plots  reveal  again  that,  for  values  of  the  ply  angle 
0  for  which  the  cross-coupling  parameter  037  reaches  a  negative  value,  smaller  displacements 
are  experienced  as  compared  to  the  case  in  which  037  is  positive. 

X  -  CONCLUSIONS 

A  dynamic  theory  of  aircraft  wings  modeled  as  thin-walled  composite  beams  of  closed 
cross-sectional  contour  was  presented.  The  theory  incorporates  a  number  of  features  essential 
for  a  reliable  prediction  of  the  free  vibration  and  aeroelastic  response  characteristics. 

A  specific  ply-angle  distribution,  inducing  a  bending-twist  cross-coupling,  was  con¬ 
sidered  and  its  influence  on  the  eigenvibration,  divergence  instability  and  static  aeroelastic 
response  was  investigated.  As  the  numerical  illustrations  reveal,  the  bending-twist  cross¬ 
coupling  term  037,  considered  in  conjunction  with  the  tailoring  technique,  plays  a  key  role  in 
the  enhancement  of  swept  wings  performance.  Physically,  this  coupling  stiffness  generates 
a  wash-out  effect,  alleviating  the  excessive  build-up  of  aerodynamic  loads.  This  is  achieved 
by  using  ply  angles  so  as  to  mininnize  or  completely  nullify  the  wash-in  effect,  as  well  as  the 
detrimental  transverse  shear  flexibility  effect.  A  similar  conclusion  can  be  drawn  in  connec- 
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tion  with  the  warping  restraint  effect,  which  is  always  present  and  should  be  used  to  enhance 
the  aeroelastic  characteristics  of  wing  structures. 

The  various  results  obtained  here  are  expected  to  contribute  to  the  understanding  of 
the  role  played  by  a  number  of  nonclassical  factors,  which,  as  shown,  can  affect  in  a  complex 
way  the  behavior  of  aircraft  wings  made  of  composite  materials. 
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Appendix 


The  reduced  rigidity  quantities  are  as  follows; 

All  =  A22  - — ,  Ai2  =  A26 - ;; - =  A21,  A13  =  2Ai2  — 

All  All  ^ 

Ai4  =  522-^%^  =  A41,  A22  =  A66-4^,  A23  =  2A22^ 

All  All  ^ 

D  Ai65i2  T'  Ty'  r,r'  /  a 

A24  =  -826 - ;; - =  A42,  A43  =  2A24-3-  (A.l) 

^11  P 

V  -  n  rr  -  P  -®i6Ai2  _  d  ^leAie 

ii44  —  JJ22 - Z — ,  2451  —  ±J26 - Z - 1  -^52  —  -D66 - Z - 

All  All  All 

Ks3=2Ks2^,  A:s4  =  fl26-%^ 

/5  All 

where  A,y,  Bij  and  Dij  denote  local  stretching,  stretching-bending  coupling  and  bending 
rigidity  quantities,  respectively. 

The  inertia  terms  are 


h  =  ^{uo  -  y  '^)mods,  h  =  ^{vo  +  x^)mods 

h  =  ^ (wo  +  x6y  +  yOx  -  Fu'^')mods,  h  =  ^ +  y^)"<i>  -  yuo  +  xvo]mods 
/s  =  ^[xibo  +  x'^dy  +  xyOx  —  xFuj{s)'^']mods 
+  j>  [i^'Oy  +  m2ds 

h  =  j>  [ywQ  +  y^Bx  +  ^yBy  -  yA^(s)^']mods 

-|-  ^  4-  m£6y  —  ma(s)^^j  m2ds 

h  =  ^ [-Aa,(5)u;o  -  xFuBy  -  yFij{s)'9x  -h  Fi^{s)^'^']mQds 
+  y  1^— ma(s)^i  —  la{s)6y  4-  a^(s)^'j  m2ds 


where 


(mo, m2)  =J2L 


(A.2) 


(A.3) 


denote  the  mass  terms. 


The  rigidity  terms  are  symmetric,  aij  =  Oj,-,  and  the  terms  of  interests  in  this  paper 
have  the  expressions 


033  =  ^  (Aiij/^  4-  2yKum  4-  A44m^)  ds 
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(A4) 


037  =  ^  {yKu  4-  K^^m)  ds,  055  =  ^  {K22^^  +  ^44771^)  ds 
056  ~  ~  ^  {FuiK2\^  +  A'^240^)  ds 

066  =  ^  +  2Ki4^Fi^a  +  ds,  077  =  ^  2—^K2zds 

The  reduced  mass  terms  are  as  follows: 

(61,  64,  65,  h\Q)  =  j>mo{l,  y^,  x^,  Fl)ds 
(^14,  ^is)  =  /  m2  (m^,  ds 


(A5) 
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Abstract 

This  paper  presents  an  integrated  approach  combining  structural  tailoring  with  the  con¬ 
verse  piezoelectric  effect  for  the  purpose  of  actively  controlling  the  vibration  and  static 
aeroelastic  characteristics  of  advanced  aircraft  wings.  The  structural  model  incorporating 
a  number  of  nonclassical  features  consists  of  a  thin/thick- walled  closed  cross  section  can¬ 
tilevered  beam  whose  constituent  layers  exhibit  elastic  anisotropic  properties.  In  addition, 
a  system  of  piezoelectric  actuators  bonded  to,  or  embedded  into  the  structure  generates  a 
localized  strain  field  in  response  to  an  injected  electric  current,  thus  producing  a  change  in 
the  static  and  dynamic  characteristics  of  the  structure.  Results  reveal  that  a  combination 
of  both  techniques  can  play  a  major  role  in  enhancing  the  vibrational  and  static  aeroelastic 
response  characteristics  of  aircraft  wings. 

1.  Introduction 

Due  to  their  outstanding  properties,  such  as  high  strength/stiffness  to  weight  ratio, 
fiber- reinforced  laminated  thick/thin- walled  structures  are  likely  to  play  an  increasing  role 
in  the  design  of  advanced  aircraft  wings.  In  addition,  a  number  of  elastic  couplings  resulting 
from  anisotropy  and  ply-angle  sequence  of  the  composite  materials  can  be  exploited  so  as  to 
enhance  the  response  characteristics.  In  this  regard,  within  the  last  two  decades,  a  technique 
referred  to  as  structural  tailoring  has  been  used  with  spectacular  results.^  It  should  be  noted. 

The  authors  are  listed  in  alphabetical  order. 

*  Professor. 

**  University  Distinguished  Professor.  Fellow  AIAA. 

***  Post-Doctoral  Research  Associate. 


1 


however,  that  structural  tailoring  is  a  passive  design  technique,  in  the  sense  that  the  control 
law  is  fixed  in  terms  of  the  considered  constitutive  equations.  This  implies  that  the  structure 
cannot  respond  adaptively  to  changes  in  its  parameters  or  external  stimuli.  To  overcome  this 
shortcoming,  additional  capabilities  must  be  built  into  the  structure.  This  is  particularly 
true  in  view  of  the  fact  that  future  generations  of  flight  vehicles  are  likely  to  operate  under 
increasingly  severe  conditions. 

An  approach  showing  good  promise  is  adaptive  control,  which  amounts  to  the  integra¬ 
tion  of  adaptive  materials  possessing  sensing  and  actuating  capabilities  into  the  structure.^"* 
Piezoelectric  materials  are  excellent  candidates  for  the  role  of  sensors  and  actuators.  In  con¬ 
trast  to  passive  structures,  in  which  the  vibrational  and  aeroelastic  response  characteristics, 
are  predetermined,  in  adaptive  structures  these  characteristics  can  be  altered  in  a  known  and 
predictable  manner.  These  adaptive  capabilities  can  be  used  to  prevent  structural  resonance 
and/or  any  other  type  of  instability,  as  well  as  to  improve  the  static  and  dynamic  response 
of  the  structure. 

In  this  paper,  the  task  of  improving  the  static  aeroelastic  response  and  free  vibration  char¬ 
acteristics  of  aircraft  wings  made  of  advanced  composite  materials  is  accomplished  through 
the  synergistic  effect  of  combining  structural  tailoring  and  adaptive  control  techniques.  The 
structural  wing  model  consists  of  a  thin/ thick- walled  closed  cross-sectional  cantilevered  beam 
whose  constituent  layers  feature  elastic  anisotropic  properties.  The  adaptive  control  capa¬ 
bility  is  achieved  by  a  voltage  feedback  via  the  converse  piezoelectric  effect.  The  induced 
localized  strain  field  produces  a  change  in  the  dynamic  characteristics  of  the  structure. 

The  global  constitutive  equations  of  wing  structures  made  of  advanced  composite  ma¬ 
terials  and  incorporating  adaptive  capabilities  are  first  derived.  These  adaptive  capabilities 
are  provided  by  piezoelectric  layers  bonded  or  embedded  into  the  structure  and  serving  both 
as  sensors  and  actuators.  Then,  based  on  related  work,®  the  equations  of  motion  and  the 
associated  boundary  conditions  of  composite  adaptive  structures  are  derived.  Feedback  con¬ 
trol  laws  relating  the  applied  electric  field  to  the  mechanical  characteristics  of  the  vibrating 
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structure  are  implemented,  thus  altering  the  frequencies  and  mode  shapes  of  the  system. 
The  obtained  results  underhne  the  fact  that  simultaneous  implementation  of  tailoring  and 
adaptive  materials  technology  can  enhance  the  static  aeroelastic  response  and  dynamic  char¬ 
acteristics  of  flight  vehicle  structures  significantly. 

2.  The  Structural  Model 

A  structural  model  consisting  of  a  thin-walled  beam  of  arbitrary  cross-section  aiming  at 
simulating  the  lifting  surface  of  advanced  flight  vehicles  is  used.  Two  systems  of  coordinates, 
namely  s,z,n  and  x,  y,  2:,  are  used  to  describe  the  kinematics  of  thin- walled  beams,  as  shown 
in  Fig.  1.  The  theory  of  thin- walled  beams  used  herein  incorporates  the  following  nonclassical 
features:  i)  Anisotropy  of  constituent  material  layers,  ii)  Transverse  shear  flexibility  and  iii) 
Primary  and  secondary  warping  effects.  In  the  light  of  ii),  the  structural  model  applies  not 
only  to  thin-walled  beams,  but  also  to  thick-walled  beams.  The  theory  is  based  also  on  the 
in-plane  cross-section  nondeformability  assumption.®'® 

Consistent  with  the  above  statements,  the  components  of  the  displacement  vector  are 


expressed  as 

u{x,y,z,t)  =  uo{z,t)  -  y(f){z,t)  (la) 

v{x,y,z,t)  =  vo{z,t)  +  x(l){z,t)  [lb) 

w  (x,  y,  z,  t)  =  wo  (z,  t)  +  Ox  (z,  t)  [y  (s)  +  nm] 

+  Oy  (z,  t)  [x(s)  -1-  ni]  -  <!>'  (z,  t)  [F,i;(s)  +  na(s)]  (Ic) 

where  I  =  cos  (n,x)  and  m  =  cos  (n,y)  denote  the  direction  cosines.  In  addition 

Ox{z,t)=lyz{z,t)-VQ{z,t)  (2a) 

0j,(z,t)  =  7xz(z,t)-Uo(z,t)  (26) 


where  6^  {z,t)  and  6y  {z,t)  denote  the  rotations  about  axes  x  and  y,  respectively,  while  7j,^ 
and  7iz  denote  the  transverse  shear  strains  in  the  planes  yz  and  xz,  respectively.  Based  on 
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Eqs.  (1),  the  axial  strain  component  reduces  to 

^ZziP'J  ^5^)0  ~  *^22(^5  0  zzi.^  J 

(3) 

where 

Szz{s,z,t)  =  w'„{z,t)  +  0j,(z,t)a:(s)  +  9'j.{z,t)y{s)  -  (f>''{z,t)F^{s) 

(4a) 

and 

Szz{s,  Z,  t)  =  t)£  +  e'^{z,  t)m  -  (f>"{z,  t)a{s) 

(46) 

are  the  axial  strains  associated  with  the  primary  and  secondary  warping,  respectively,  in 

which  primes  denote  derivatives  with  respect  to  z.  The  membrane  and  transverse  shear 

strain  components  can  be  expressed  in  the  form 

Ac 

Sazis,  Z,  t)  =  Ssz{s,  Z,  t)  +  2—(j>{z,  t) 

(5) 

and 

Snz{s,z,t)  =  [ey{zj)  +  u'„{z,t)\£+  [9z{z,t)  +  v'„{z,t)]  m 

(6a) 

respectively,  where 

Saz{s,Z,t)  =  -  \9y{z,t)  +  u„{z,t)\  m  +  [9z{z,t)  -f  v'„{z,t)]  £ 

(66) 

In  the  above  equations,  Uo{z,t),  Vo{z,t)  and  Wo{z,t)  represent  rigid-body  translations 

in  the 

x,y  and  z  directions  and  4>{z,t)  represents  the  twist  about  the  z-axis,  (see  Fig.  2).  Moreover, 
h  =  h{s)  denotes  the  wall  thickness,  allowed  to  vary  in  the  circumferential  direction,  Ac  the 
cross-sectional  area  bounded  by  the  contour  midline,  j3  the  total  length  of  the  contour  midline 

and  t  the  time.  In  addition,  F^^{s)  denotes  the  warping  function  defined  as®'® 

■ 

Fu,{s)=  f  [r„(s)-'0]ds 

Jo 

(7) 

where 

(8) 
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is  the  torsional  function,  in  which  ^{•)ds  denotes  the  integral  along  the  closed  midline  contour, 
and  r„(s)  and  a  are  geometric  quantities  defined  as 

r„(s)  =  x(s)£  +  t/(s)m,  a  =  —xm  +  y£  (9) 

When  the  structure  exhibits  infinite  rigidities  in  transverse  shear,  6y{z,t)  — »■  —u'g{z,t)  and 
Ox{z,t) 

3.  Constitutive  Equations  for  Adaptive  Wings 

As  recently  conjectured, the  adaptive  capabilities  of  piezoelectric  devices  can  result 
in  superior  control  compared  to  standard  passive  or  even  active  control.  Piezoelectric  mate¬ 
rials,  integrated  into  the  load-bearing  structure  by  means  of  surface  bonding  or  embedding, 
serve  as  networks  of  actuator /sensor  systems. 

The  linear  three-dimensional  piezoelectric  constitutive  equations  can  be  expressed  in 
contracted  index  notation^^ 


£  s 

(Ti  —  CfjSj  6ki  £kt  Dk  ~  ^kjSj  "b  ^kt 


(10a,  h) 


where  a,-  and  Sj{i,j  =  1,2, ...  ,6)  denote  the  stress  and  strain  components,  respectively,  in 
which 

c._/  Spr  ioT  p  =  rj  =  1,2,3  ,  . 

^  1  25'pr  for  p  ^  r,j  =  4,5,6 

Moreover,  Cfj,  e^i  and  are  the  elastic  (measured  under  constant  electric  field),  piezoelec¬ 
tric  and  dielectric  constants  (measured  under  constant  strain),  while  and  Dk(k  =  1,2,3) 
denote  the  electric  field  intensity  and  electric  displacement  vector,  respectively.  Summation 
over  repeated  indices  is  implied  in  Eqs.  (10).  Equations  (10a)  and  (106)  describe  the  converse 
and  direct  piezoelectric  effects,  respectively. 

In  piezoelectric  adaptive  structures,  the  direct  effect  is  used  for  sensing  and  the  converse 
effect  is  used  to.  generate  control  forces.  Equations  (10)  are  valid  for  the  most  general 
anisotropic  case,  i.e.,  for  triclihic  crystals.  In  the  following,  the  piezoelectric  anisotropy  is 
restricted  to  the  case  of  hexagonal  symmetry,  the  n-axis  being  an  axis  of  rotary  symmetry 
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coinciding  with  the  direction  of  polarization^^  (thickness  polarization).  In  this  case,  the 
piezoelectric  continuum  is  characterized  by  a  reduced  number  of  elastic  piezoelectric  and 
dielectric  coefficients.  We  also  assume  that  the  electric  field  vector  £i  is  defined  by  the 
component  oiily?  E\  '=  €2  —  and  that  the  electric  field  intensity  E^  is  uniform  over 
the  entire  area.  Due  to  the  fact  that  the  voltage  is  distributed  uniformly,  we  assume  that 
in  the  static  case  E^  is  constant  while  in  the  dynamic  case  it  depends  on  time  alone.  The 
piezoelectric  actuators  are  distributed  over  the  entire  span  and  their  distribution  in  the 
circumferential  and  transversal  directions  is  defined  by 


Rk{n)  =  H[n-  -  H  [n  -  nf^k+))  (12a) 

Rkis)  =  H(^s-  -H  (^s-  (126) 


in  which  H  denotes  the  Heaviside  unit  step  function.  It  is  assumed  that  the  host  structure 
and  the  actuator  consist  of  t-  layers  exhibiting  general-type  orthotropy,  and,  respectively  of 
p  piezoelectric  layers  exhibiting  hexagonal  symmetry. 

Before  we  derive  the  global  one-dimensional  constitutive  equations  for  adaptive  wing 
structures,  it  will  prove  convenient  to  derive  their  two-dimensional  counterparts.  Integrating 
the  actual  three-dimensional  constitutive  equations  over  the  wall  thickness  and  assuming 
that  the  hoop  stress-resultant  is  negligibly  small  when  compared  with  the  remaining 
stresses,  the  two-dimensional  constitutive  equations  are 


N zn  —  ^44452 


(14) 

(15) 
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where  Nzz  and  denote  tangential  stress  resultants,  Ngn  denotes  the  transverse  shear 
stress  resultant  and  Lzz  and  Lsz  denote  the  stress  couples,  all  quantities  depending  on  s,z 
and  t.  Moreover,  Kij  denote  the  modified  local  stiffness  coefficients,  listed  in  the  Appendix, 
and  N^z  and  L^z  denote  the  piezoelectrically  induced-stress  resultant  and  stress  couple, 
respectively,  expressed  as 

^z“z(^>0=  (l-^)  («(ib+)  -  «(ib-))  4\^R{k)is)  (16a) 

Llzis,  t)=Yl  4^\i)R{k){s)  (n(jt+)  -  n(i-))  ^ (n(jt+)  +  n(fc-))  -  ^  (166) 

Jt=i 

where,  in  the  case  of  symmetrically  located  piezoactuators,  the  term  underscored  vanishes. 

4.  The  Equations  of  Motion  for  Adaptive  Wing  Structures 

From  Ref.  5,  the  one-dimensional  stress-resultants  and  stress  couples  are  given  by 

Tz{z,t)  =  ^  Nzzds,  Qx{z,t)  =  j>  {-Nsz'm  + Nzn^)ds  (17a,  6) 

Qy{z,t)  {N3zi  +  Nznm)ds,  Mx{z,t)  =  j>  {yNzz  ■{■LzzTn)ds  (17c,  d) 

My{z,t)=j{xNzz-\-Lzzf)ds,  Mz{z,t)  =  ^  Nsz^^ds  (17e,/) 

and 

Bu,{z,t)  =  j)  [Fuis)Nzz  +  a{s)Lzz]ds  {17g) 

where  Tz  is  the  axial  force,  Qx  and  Qy  are  shear  forces,  Mx,  My  and  Mz  are  the  bending  and 
twist  moments,  respectively,  while  is  a  bimoment  global  quantity.  In  view  of  Eqs.  (13) 
and  (15),  the  stress  quantities  Tz,  Mx,  My  and  B^^  can  be  cast  in  the  more  convenient  form 

Tz  =  fz-  fz,  Mx  =  Mx-  Mx,  My  =  My-  ky,  B^  =  B^-  K  (18a-d) 

where  single  and  double  overcarets  identify  purely  mechanical  and  piezoelectrically  induced 
terms,  respectively.  The  latter  terms  are 

Tz  =  ^  X)  4^^^(it)(^)  ("(it+)  -  "(fc-))  4^1  [\  {^(k+)  +  "(it-)) 
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It  is  readily  seen  from  Eqs.  (19a-d)  that  the  piezoelectrically  induced  stress  resultants  are 
proportional  to  the  injected  electric  current  £z-  In  the  case  of  actuators  placed  symmetrically 
throughout  the  thickness  of  the  beam,  the  underlined  terms  in  Eqs.  (19)  vanish. 

One  reason  for  employing  advanced  composite  materials  in  flight  vehicle  design  lies  in 
the  fact  that  they  permit  the  use  of  specific  lay-up  and  fiber  orientations  so  as  to  induce 
preferred  elastic  couplings  having  enhanced  effects  on  the  response  characteristics.  As  studies 
on  vibration  and  subcritical/critical  aeroelastic  behavior  of  wing  structures  reveal, 
bending-twist  coupling  plays  a  major  role.  Additional  beneficial  effects  of  this  cross  coupling 
have  been  discussed  recently  in  Ref.  20.  In  the  case  of  wing  structures  modeled  as  thin-walled 
beams,  the  ply-configuration  inducing  such  an  elastic  coupling  is  referred  to  in  Refs.  12  as 
the  circumferentially  antisymmetric  stiffness  configuration  and  in  Ref.  13  as  the  symmetric 
configuration.  The  associated  ply-angle  distribution  is  governed  by  the  law  (see  Fig.  3) 

e{y)  =  -e{-y)  (20) 

The  resulting  equations  for  adaptive  wing  structures  characterized  by  the  above  men¬ 
tioned  ply-angle  configuration  are: 

-  ae64>""  +  073^1  —  a65(i’o'  +  O'x)  +  +  ruz  =  (64  -f  65)  ^  -  {pig  ffpig)3"  (21a) 
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(216) 


055  (v'o'  +  ^x)  +  056^"^  + Py  =  biVo 

033^x  +  037^”  ~  (^°  ~ 

as  well  as  by  the  boundary  conditions  at  the  root  (z  —  0) 

<^  =  0,  =  0  ,  uo  =  0,  01  =  0  (22a-d) 

and  at  the  tip  (z  =  L) 

-  a66<f>"'  +  07301  +  077^6'  =  -  (610  +  613)  5  a66<A”  +  “65  (Uq  +  0x)  =0  (23a,  6) 

055  (^0  +  0i)  +  a5s(f>"  =  0,  0330^  +  az7<j>  =  Mx  (23c,  d) 

The  coefficients  appearing  in  Eqs.  (21)  and  (23)  stand  for  the  stiffness  and  mass  quantities 
and  are  given  in  the  Appendix.  The  existence  of  two  bending-twist  stiffness  coupling  terms, 
namely,  037  =  073  and  ase  =  055,  should  be  noted.  The  latter  is  induced  by  the  warping 
restraint  effect  and  its  influence  is  in  general  somewhat  less  important  than  that  of  the 
former.  Their  variation  with  the  ply-angle  is  presented  in  Fig.  16  of  Ref.  5.  The  singly  and 
doubly  dashed  underhned  terms  in  Eqs.  (21)-(23)  are  associated  with  the  warping  restraint 
and  warping  inertia,  respectively.  Because  the  piezoactuators  are  distributed  over  the  entire 
wing  span,  derivatives  with  respect  to  z  of  piezoelectrically-induced  terms  in  Eqs.  (21)  vanish. 
This  explains  why  their  contribution  does  not  appear  in  the  equations  of  motion;  it  appears 
only  as  a  nonhomogeneous  term  in  the  boundary  conditions. 

It  should  be  observed  here  that  for  the  free  vibration  problem  the  transverse  load  py  and 
the  twisting  moment  must  be  ignored  in  Eqs.  (21)  and  for  the  static  aeroelastic  problem 
the  inertia  terms  must  be  omitted  from  Eqs.  (21)  and  (23).  Consistent  with  strip-theory 
aerodynamics,  the  unsteady  lift  force  py  and  torsional  aerodynamic  moment  rux,  both  per 
unit  span,  can  be  expressed  as 

Py{z)  =qncao  (^1^0  +  <l>-  v'o  tan  A)  -  NW/2L 
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(24a,  b) 


mt{z)  =qnCaoe  (^(j)o  +  <l>-  v'g  tan  a]  +  qnC^CMAC  -  NWdf2L 

Here  denotes  the  dynamic  pressure  normal  to  the  leading  edge  of  the  swept  wing, 

c  the  chord  of  the  wing,  Og  the  “corrected  lift”  curve  slope  coefficient,  A  the  angle  of  sweep 
(considered  positive  for  swept-back),  e  the  offset  between  the  aerodynamic  and  reference 
axis,  <^o  the  rigid  angle  of  attack  (measured  in  planes  normal  to  the  leading  edge).  Cm  AC  the 
wing  section  pitching  moment  coefficient  (whose  influence,  as  usual,  is  disregarded),  WI2L 
the  weight  per  unit  length  of  wing  and  N  the  load  factor  normal  to  the  wing  surface,  whose 
expression  is 

N  =  J  (j)o  +  <f>-Vg  tan  a)  dz  (25) 

The  static  aeroelastic  response  is  analyzed  both  in  the  subcritical  range,  i.e.,  for  velocities 
qn  <  (9n)D)  and  in  the  critical  case  as  well,  where  {qn)D  denotes  the  divergence  dynamic 
pressure.  As  a  general  remark,  Eqs.  (24)  reveal  that  for  A  <  0,  i.e.,  for  swept-forward 
wings,  the  aeroelastic  bending- twist  coupling  results  in  an  increase  in  Py(z)  and  m2(z), 
which  in  turn  results  in  a  deterioration  of  the  subcritical  static  aeroelastic  response  and 
in  a  dramatic  decrease  of  the  divergence  speed.  Whereas  the  goal  of  the  subcritical  static 
aeroelastic  analysis  consists  of  the  determination  of  the  distribution  of  the  effective  angle 
of  attack  c^eff  ^iid  of  the  lift  force,  cis  affected  by  the  elastic  deformations,  the  study  of 
the  critical  case  involves  the  determination  of  divergence  instability  conditions.  Clearly,  the 
main  target  of  tailoring  applied  to  swept-forward  wings  is  to  yield  a  decrease  in  the  effective 
angle  of  attack,  and  implicitly  in  the  aeroelastic  lift  and,  as  a  result,  an  increase  in  the 
critical  divergence  speed.  Whereas  the  study  of  the  subcritical  static  aeroelastic  response 
requires  the  solution  of  the  integral-differential  system  of  equations  obtained  by  inserting 
Eqs.  (24)  and  (25)  into  Eqs.  (21),  the  determination  of  the  divergence  speed  leads  to  the 
solution  of  an  eigenvalue  problem,  where  the  divergence  speed  plays  the  role  of  eigenvalue. 
Structural  tailoring  applied  to  the  vibration  of  wing  structures  must  result  in  an  increase 
in  the  eigenfrequencies  without  weight  penalties.  The  determination  of  natural  frequencies 
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requires  the  solution  of  an  eigenvalue  problem.^^ 

In  spite  of  the  mathematical  complexities  involved  in  the  solution  of  the  above  men¬ 
tioned  problems,  in  which  the  eigenvalue  appears  both  in  the  differential  equations  and  the 
boundary  conditions,  the  computational  methodologies  used  proved  to  be  extremely  power¬ 
ful.  The  spatial  Laplace  transform  method  yields  exact  solutions, but  is  computationally 
laborious.  On  the  other  hand,  the  extended  Galerkin  method  yields  approximate  solutions 
in  excellent  agreement  with  the  exact  ones  and  with  less  computational  effort. 

5.  The  Control  Law  and  the  Closed-Loop  Eigenvalues 

For  feedback  control,  the  applied  electric  field  £3,  upon  which  the  piezoelectrically  in¬ 
duced  moment  depends,  is  a  function  of  the  wing  motion.  Two  simple  control  laws  are  being 
considered.  The  first  control  law,  denoted  by  CLl,  requires  that  the  injected  electric  field 
Sz  is  proportional  to  the  bending  moment  Afx(O)  at  the  wing  root,  which  implies  that 

^3  =  (26) 

Upon  considering  the  boundary  condition  (23d),  as  well  as  Eq.  (19b),  the  control  law  given 
by  Eq.  (26)  can  be  rewritten  as 

«i(L)+/3<!S'(i)-tX(0)  =  0  (27) 

where  fz  =  037/033  and  kp  denotes  the  feedback  gain.  This  control  law  expresses  the  fact  that 
the  bending  moment  at  the  wing  tip  induced  by  piezoelectric  strain  actuation  is  proportional 
to  the  mechanical  bending  moment  at  the  wing  root.  This  control  law  was  used  in  Refs.  14- 
16.  The  second  control  law,  denoted  by  CL2,  requires  that  the  applied  electric  field  S3  is 
proportional  to  the  vertical  deflection  of  the  beam  tip  vo[L).  This  results  in  the  condition 

e',{L)  +  f3cf>'{L)  =  kpVo{L)  (28) 

Equation  (28)  expresses  the  fact  that  the  boundary  moment  control  at  the  wing  tip,  induced 
by  piezoelectric  strain  actuation,  is  proportional  to  the  transversal  deflection  at  the  wing 
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tip.  Control  law  (28)  is  similar  to  the  one  proposed  in  Ref.  24.  It  should  be  noticed  that  the 
feedback  ^ains  kp  and  kp  are  nondimensional  and  dimensional  quantities,  respectively.  The 
non  dimensional  counterpart  of  kp  is  kp  =  kpL^. 

6.  Numerical  Examples 

The  equations  derived  here  are  general,  in  the  sense  that  they  are  valid  for  a  beam  of 
arbitrary  cross  section,  as  well  as  for  piezoelectric  actuators  arbitrarily  located  throughout 
the  waU  thickness  and  along  the  circumference  of  the  beam.  However,  in  the  present  case 
a  biconvex  profile,  typical  of  supersonic  wing  airplanes,  is  adopted.  It  is  assumed  that 
the  piezoceramic  actuators  used  here  are  mounted  symmetrically  on  the  upper  and  bottom 
surfaces  of  the  wing.  The  host  structure  is  assumed  to  be  of  a  graphite/epoxy  material  whose 
elastic  characteristics  are 

El  =  30  X 10®  psi,  Et  =  0.75  x  10®  psi 
Git  =  0-37  X 10®  psi,  Gtt  =  0-45  x  10®  psi 
HTT  =  y-LT  =  0.25,  p  =  14.3  X  10~®  Ibsec^/in'* 

where  the  subscripts  L  and  T  denote  directions  parallel  and  transversal  to  the  fibers,  re¬ 
spectively.  The  geometrical  wing  characteristics  are  displayed  in  Fig.  1  and  the  piezoelectric 
actuators  are  made  of  PZT-4  ceramic,  whose  properties  are  given  in  Ref.  17. 

The  associated  differential  eigenvalue  problem  has  been  discretized  in  space  by  the  ex¬ 
tended  Galerkin  method.^^  The  obtained  results  were  checked  by  means  of  an  exact  solution 
based  on  a  spatial  Laplace  transform^^’^^  and  the  agreement  was  excellent. 

Tailoring  and  adaptive  control  are  applied  simultaneously  to  enhance  the  vibrational 
characteristics  of  wing  structures.  Results  are  displayed  in  Figs.  4  and  5  and  Figs.  7  and 
8  for  CLl  and  CL2,  respectively.  It  should  be  mentioned  that  for  0  /  0  and  0  ^  90° 
the  bending-twist  coupling  causes  the  adaptive  control  to  be  somewhat  weaker  in  places 
where  the  twist  is  very  strong.  On  the  other  hand,  for  9  =  0°  and  90°,  where  the  bending 
decouples,  the  adaptive  control  becomes  very  effective.  As  a  general  rule,  the  frequencies 
were  normalized  with  respect  to  the  ones  corresponding  to  the  uncontrolled  case  kp  =  0,  and 
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to  on-axes  ply  angle  configuration,  9  =  0.  The  trend  in  the  eigenfrequencies  revealed  in  this 
paper  using  CL2  agrees  fully  with  one  obtained  theoretically  and  verified  experimentally.^^ 
As  can  be  concluded  from  Figs.  5  and  8,  the  first  eigenfrequency  changes  gradually  with 
the  ply  angle  from  one  corresponding  to  a  pure  bending  mode  (labelled  as  5)  at  0  =  0  to 
a  coupled  bending-twist  mode  for  0  <  ^  <  90°  to  a  pure  bending  mode  for  6  =  90°.  On 
the  other  hand,  the  second  eigenfrequency  changes  from  a  pure  bending  mode  at  0  =  0  to 
a  pure  twist  mode  (labelled  as  T)  for  6  =  90°,  and  the  third  eigenfrequency  changes  from 
pure  twist  for  0  =  0  to  pure  bending  for  6  =  90°. 

Plots  of  the  first  closed-loop  mode  are  shown  in  Fig.  6  for  three  values  of  feedback  gains. 
It  is  easy  to  see  that  the  original  uncontrolled  modes  differ  considerably  from  the  controlled 
modes.  In  Figs.  9a  and  9b,  the  reduced  divergence  speed  parameter  {qnCao)j)  is  displayed  as  a 
function  of  6  for  various  values  of  the  feedback  gain  using  CLl  and  for  the  cases  of  a  straight 
wing,  A  =  0,  and  a  swept-forward  wing,  A  =  —30°.  As  readily  seen,  the  use  of  both  tailoring 
and  adaptive  control  with  negative  feedback  gains  result  in  a  significant  increase  in  the 
divergence  speed.  In  Figs.  10a,  10b  and  11,  the  normalized  effective  angle  of  attack 
is  displayed  for  nonadaptive  wings,  =  0,  and  for  adaptive  wings,  <  0,  respectively  . 
The  results  reveal  that  a  strong  attenuation  of  aeroelastic  loads  is  experienced  for  kp  <  0,  as 
compared  to -the  case  corresponding  to  kp  =  0.  As  Fig.  11  reveals,  as  the  magnitude  of  the 
feedback  gain  increases,  the  attenuation  becomes  more  and  more  pronounced.  It  should  be 
noted  that  in  the  case  of  CL2,  improved  vibrational  characteristics  are  obtained  for  positive 
feedback  control  gains,  kp  >  0. 

The  results  demonstrate  that  the  use  of  both  tailoring  and  adaptive  control  techniques 
can  yield  dramatic  improvements  in  the  dynamic  and  static  aeroelastic  characteristics  of 
aircraft  wings. 

7.  Summary  and  Conclusions 

A  comprehensive  structural  model  of  aircraft  wing  modeled  as  a  thin/thick- walled  beam 
of  arbitrary  closed  cross  section,  cantilevered  aircraft  wing  structures  incorporating  adaptive 
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capabilities  via  the  inverse  piezoelectric  effect  was  developed.  The  anisotropic  wing  model 
incorporates  a  number  of  nonclassical  features  characterizing  actual  wing  structures. 

The  unique  directional  properties  of  composite  materials  permits  tailoring,  which 
amounts  to  enhancing  wing  structure  response,  whereeis  incorporation  of  piezoactuators  per¬ 
mits  vibration  and  static  aeroelastic  control  of  the  structure.  It  is  demonstrated  here  that 
the  synergistic  combination  of  tailoring  of  anisotropic  composite  materials  and  control  by 
means  of  adaptive  materials  results  in  a  wing  with  better  dynamic  and  static  aeroelastic 
characteristics  than  would  have  resulted  from  tailoring  or  control  alone.  Incorporation  of 
both  technologies  can  provide  an  expanded  performance  envelope  of  flight  vehicles,  without 
weight  penalties. 
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Appendix 


The-modified  local  stiffness  coefficients  are 

Ku=A-a-^,  Ku  =  A2(-^^  =  Kiu  Kn=2Kuj  (Al-3) 

Ku  =  Bi2-^^^  =  K»,  K22  =  Ass  -  K23  =  2K22^  (A4.6) 

All  All  P 

K2a  =  B2e-^^^  =  K42,  K^,=2K24^,  K44  =  D22-^  (A7-9) 

All  P  All 

Ka  =  K32  =  B6S-%^,  K33  =  2K32^,  Km  =  B2S-%^  (AlO-13) 

All  All  P  All 

where  A,j,  Bi,  and  Dij  denote,  respectively,  the  stretching,  bending-stretching  and  bending 
stiffness  quantities  associated  with  the  entire  structure,  where  the  structure  consists  of  TV  = 
%  +  p  host  and  piezoactuator  layers. 

The  stiffness  and  mass  coefficients  have  the  form 
an  =  /  Kiids,  ai4  =  041  =  —  ji  TV^mds,  022  =  ^  [Kii  +  2xK\4l  -f  7^44^^)  ds  (A14-16) 

033  =  ^  {Kny'^  +  2yKi4m  +  K44m?^  ds,  037  =  073  =  ^  {yl<i3  +  I<43rn)  ds  (A17, 18) 

044  =  y  (JK221^^  +  A447^^  ds,  055  =  ^  {K22^^  +  A44m^ j  ds  (A19, 20) 

^56  =  065  ~  ~  ^  {FuK2li  +  7^2407)  ds  (-^21) 

066  =  ^  (KiiFu  +  27Vi4Fi^a  +  7^440^)  ds,  ajj  =  ^  2-^K2zds  (A22, 23) 

(61,64,65,610)  =  y  mo  (1,  j/^,  x^,  FI,')  ds,  (614, fei5, tis)  =^m2  (m^, 7^, a^)ds  (A24,25) 

where 

(mo,  m2)  =  P(k)  dn  (A26) 
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FIGURE  CAPTIONS 


Fig.  T-  The  Thin- Walled  Cantilever  Beam 

Fig.  2-  Displacement  Field  and  Coordinate  Systems 

Fig.  3-  Circumferentially  Antisymmetric  Stiffness  Configuration,  6{y)  =  —0{—y) 

Fig.  4-  First  Normalized  Eigenfrequency  (iJi  =  uiijufj)  Versus  the  Feedback  Gain  Using 
CLl  for  Various  Ply  Angles 

(wjv  =  44  rad/s,  corresponding  to  A:p  =  0  and  0  =  0) 

Fig.  5-  Three  Normalized  Bending-Twist  Coupled  Eigenfrequencies  (cU,-  =a;,/a;jv)  Versus 
the  Ply  Angle  for  No  Control  and  One  Feedback  Gain  Using  CLl 
{uff  =  44  rad/s,  corresponding  io  kp  =  Q  and  0  =  0) 

Fig.  6-  First  Normalized  Mode  for  Three  Feedback  Gains  Using  CLl-,  0  =  0 

Fig.  7-  First  Normalized  Coupled  Eigenfrequency  (aJi  =  versus  Feedback  Gain 

Using  CL2  for  Several  Values  of  0;  =  44  rad/s 

Fig.  8-  Three  Normalized  Coupled  Eigenfrequencies  (aJ,  =  uiifiON)  for  Various  Ply  Angles 
and  No  Control  and  One  Feedback  Gain  Using  CL2-,  Wjv  =  44  rad/s 

Fig.  9a-  The  Divergence  Speed  Parameter  Versus  the  Ply  Angle  for  No  Control  and  One 
Feedback  Gain  Using  CLl]  A  =  0 

Fig.  9b-  The  Divergence  Speed  Parameter  Versus  the  Ply  Angle  for  No  Control  and  One 
Feedback  Gain  Using  CLl;  A  =  —30° 

Fig.  10a-  The  Normalized  Effective  Angle  of  Attack  as  a  Function  of  the  Position  Along  the 
Wing  Span  for  Three  Sweep  Angles;  fcp  =  0,  0  =  45°,  135°  (^flight  =  2psi) 

Fig.  10b-  The  Normalized  Effective  Angle  of  Attack  as  a  Function  of  the  Position 
Along  the  Wing  Span  for  Three  Sweep  Angles  Using  CLl]  kp  =  0.25,  0  = 

45°,  135°  (^flight  =  2psi) 

Fig.  11-  The  Normalized  Effective  Angle  of  Attack  as  a  Function  of  the  Position  Along 
the  Wing  Span  for  Various  Feedback  Gains  Using  CLl]  A  =  —15°,  0  = 
45°  (^flight  =  2psi) 
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Fig.  3  -  Circumferentially  Antisymmetric  Stiffness  Configuration,  9{y)  =  —9{—y) 


Fig.  4-  First  Normalized  Eigenfrequency  (oJi  =a;i/u;^)  Versus  the  Feedback  Gain  Using 
CLl  for  Various  Ply  Angles 

=  44  rad/s,  corresponding  to  fcj,  =  0  and  0  =  0) 


9  (degree) 


Fig.  5-  Three  Normalized  Bending-Twist  Coupled  Eigenfrequencies  (u7,-  =  Versus 

the  Ply  Angle  for  No  Control  and  One  Feedback  Gain  Using  CLl 
=44  rad/s,  corresponding  to  /cp  =  0  and  0  =  0) 


Fig.  6-  First  Normalized  Mode  for  Three  Feedback  Gains  Using  CTl;  0  =  0 


Fig.  7-  First  Normalized  Coupled  Eigenfrequency  (wi  =wi/u;;\r)  versus  Feedback  Gain 
Using  CL2  for  Several  Values  of  9;  wjy  =  44  rad/s 


Ply  angle  (degree) 


Fig.  8-  Three  Normalized  Coupled  Eigenfrequencies  (ui  =  for  Various  Ply  Angles 

and  No  Control  and  One  Feedback  Gain  Using  CL2\  =  44  rad/s 


Fig.  9a-  The  Divergence  Speed  Parameter  Versus  the  Ply  Angle  for  No  Control  and  One 
Feedback  Gain  Using  CL\\  A  =  0 


Fig.  9b-  The  Divergence  Speed  Parameter  Versus  the  Ply  Angle  for  No  Control  and  One 
Feedback  Gain  Using  CTl;  A  =  —30° 


Fig.  10a-  The  Normalized  Effective  Angle  of  Attack  as  a  Function  of  the  Position  Along  the 
Wing  Span  for  Three  Sweep  Angies;  kp  =  0,  6  =  45°,  135°  (^flight  =  2psi^ 


Fig.  10b-  The  Normalized  Effective  Angle  of  Attack  as  a  Function  of  the  Position 
Along  the  Wing  Span  for  Three  Sweep  Angles  Using  CLl]  kp  =  0.25,  0  = 

45°,  135°  (^flight  =  2psi) 


Fig.  11-  The  Normalized  Effective  Angle  of  Attack  as  a  Function  of  the  Position  Along 
the  Wing  Span  for  Various  Feedback  Gains  Using  CL\\  A  =  -15°,  6  = 
45°  (^flight  =  2psi) 


STRUCTURAL  MODELING  OF  LOW-ASPECT  RATIO  COMPOSITE  WINGS* 
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SUMMARY 

This  paper  is  concerned  with  the  aeroelastic  tailoring  of  a  structural  model  consisting 
of  a  rigid  fuselage  and  a  low-aspect  ratio  wing  made  of  composite  materials.  The  wing  is 
modeled  as  a  trapezoidal  plate  with  root  and  tip  chords  paraUel  to  the  flow  and  with  general 
sweep.  The  fuselage  is  capable  of  plunge  and  pitch  and  the  elastic  wing  model  includes  shear 

deformations  but  ignores  rotatory  inertia. 

1.  INTRODUCTION 

Recently,  and  particularly  in  the  past  decade,  solutions  to  the  flutter  problem  based  on 
anisotropic  wings  made  of  composite  materials  have  been  presented  by  many  authors.  Much 
of  this  work  refined  the  concept  of  “composite  tailoring”  first  put  forth  by  Krone^  During 
the  same  time  frame,  many  researchers  have  come  to  the  conclusion  that  the  divergence 
problem  associated  with  forward  swept  wings  is  in  fact  a  body-freedom  flutter  phenomenon. 
Body-freedom  flutter  is  the  coupling  of  a  flexible  aircraft  mode  with  a  rigid-body  mode 
due  to  interaction  with  aerodynamic  loads.  One  common  example  of  this,  and  the  one 
receiving  attention  here,  is  the  coupling  of  flexible  wing  bending  and  rigid-body  pitch.  Work 
by  Weisshaar  et  al.^  leads  to  the  conclusion  that  a  clamped  wing  model  generally  predicts 
larger  improvements  in  the  instabiUty  speed  for  composite  tailoring  than  for  models  that 
include  rigid  body  motions.  An  accurate  prediction  of  the  critical  speed  in  the  presence 
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of  structural  paxameter  vaxiations  must  be  a  part  of  any  structural  model  aimiug  to  model 
general  planforms  accurately. 

The  requirements  of  modern  wing  planforms  place  additional  constraints  on  structural 

models.  In  particular,  many  wings  of  practical  interest  are  of  too  low  an  aspect  ratio  to  be 

approximated  by  a  beam  model.  Similarly,  the  sweep  angle  should  not  be  confined  to  aft 

only,  particularly  for  an  optimization  problem.  Two-dimensional  models  capable  of  accepting 

fiber  directional  information  is  suggested  in  this  case,  but  examples  of  such  models  in  the 

•  •34 

literature  concerned  with  this  application  are  rare  ’  . 

The  structural  model  suggested  by  the  above  physical  requirements  is  a  two-dimensional 
anisotropic  plate  consisting  of  several  variable  thickness  layers  of  generally  orthotropic  ma¬ 
terial.  Because  the  interest  lies  in  composite  materials,  it  is  also  important  to  include  shear 
deformability  in  the  model.  Fuselage  rigid-body  degrees  of  freedom  must  also  be  included. 
Two  compromises  on  generality  contribute  to  the  goal  of  low  computational  effort  without 
significantly  affecting  the  desire  to  include  all  the  relevant  physical  effects.  Realistic  amounts 
of  airfoil  camber  have  little  effect  on  displacements  from  unsteady  pressure  loading,  suggest¬ 
ing  that  camber  may  be  omitted  for  aeroelastic  problems.  Similarly,  the  conclusion  that 
rotatory  inertia  terms  for  wings  of  practical  thickness  can  be  ignored  is  supported  by  the 

text  by  Librescu®. 

The  interest  lies  in  a  structural  model  that  includes  the  essential  physical  characteristics 
required  for  flutter  analysis  of  modern  low-aspect  ratio  composite  wings.  The  model  should 
be  sufficiently  accurate  to  account  for  all  important  structural  parameters  and  yet  sufficiently 
simple  so  as  not  to  be  computationaUy  intensive.  From  the  foregoing  considerations,  a  struc¬ 
tural  model  satisfying  the  imposed  requirements  has  emerged.  The  model  is  a  trapezoidal 
plate  with  root  and  tip  chords  parallel  to  the  flow.  There  are  2k  symmetrically  stacked, 
variable  thickness,  generally  orthotropic  laminae  in  the  laminate  (Figs.  1  and  2).  Mindhn 
shear  deformabihty®  is  included  with  a  shear  correction  factor  of  1.0.  The  wing  is  attached 
to  a  rigid  fuselage  capable  of  pitch  and  plunge.  This  is  the  simplest  model  retaining  the 
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essential  physical  chaxacteristics  of  low-aspect  ratio  composite  wings  with  general  sweep.  It 
is  still  an  involved  formulation,  including  three  displacement  variables  and  natural  boundary 
conditions  compUcated  by  the  sweep  angles.  To  produce  a  reasonably  accurate  solution  with 
as  few  terms  as  possible,  the  series  solution  is  in  terms  of  quasi- comparison  functions,  which 

guarantees  fast  convergence  . 

2.  DEFINITION  OF  DISPLACEMENTS  AND  BOUNDARY  CONDITIONS 

The  mathematical  model  under  consideration  consists  of  a  rigid  fuselage  attached  to  a 
flexible  wing.  The  fuselage  is  assumed  to  undergo  two  rigid-body  displacements,  namely, 
plunge  and  pitch.  On  the  other  hand,  the  wing  is  assumed  to  act  as  an  elastic  plate  clamped 
at  the  root  and  free  at  the  the  other  three  boundaries  (Fig.  1).  The  transverse  displacement 
of  a  typical  point  in  the  fuselage  has  the  form 

Wf{t)  =  luc(t)  +  (a:  -  xc)‘tpc{t) 

where  wc{t)  is  the  plunge,  defined  as  the  transverse  displacement  of  the  system  (fuselage 
and  wing)  mass  center  C  and  rpc  is  the  pitch,  defined  as  the  rotation  of  the  system  about 
an  axis  parallel  to  y  and  passing  through  C,  and  ic  is  the  distance  from  the  origin  of  xyz 
to  point  C.  Moreover,  the  displacements  of  a  typical  point  in  the  plate  are  as  follows; 


u{x,y,z,t)  =  uo(x,y,t)  +  ^^i(a^»!/)0 

(16) 

v{x,y,z,t)  =  VQ{x,y,t)  +  z‘tl)y{x,y,t) 

(Ic) 

Wp{x,y,t)  =  wc(t)  +  (x  —  xc)tpc{^)  + 

(Id) 

where  uq  and  uq  are  midplane  elastic  deflections  in  the  i  and  y  directions,  respectively,  i/>x 
and  T/.J,  are  angular  displacements  of  a  fine  normal  to  the  nominal  plane  of  the  plate  due  to 
elasticity  and  w  is  the  elastic  part  of  the  transverse  displacement. 

Because  the  wing  is  attached  to  a  fuselage  assumed  to  be  rigid,  tn  =  0,  i/’x  =  0  and 
ipy  =  0  along  y  =  0.  These  are  the  only  geometric  boundary  conditions  of  the  problem. 
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3.  THE  EXTENDED  HAMILTON’S  PRINCIPLE 

The  dynamical  problem  can  be  formtilated  by  means  of  the  extended  Hamilton’s  princi¬ 
ple,  which  can  be  written  in  the  form® 

(^ST  -SV  +  Wnc)  dt  =  0,  Swc  =  =  Sw  =  =  5i/>y  =  0  at  t  =  tut2  (2) 

where  T  is  the  kinetic  energy,  V  is  the  potential  energy  and  8W nc  is  the  virtual  work 
performed  by  the  nonconservative  forces. 

The  kinetic  energy  arises  from  two  sources,  the  motion  of  fuselage  and  of  the  wing,  and 
has  the  expression 

T  =Tf  -H  Tu,  =  ^  y  [ihc  +  (a:  -  ^c)  ^c]  drrif  +  2  Jq  mwjdxdy 
=  -  [Mc^c  +  ^ci’c)  +  mwdxdy 

^  lo  Jle  ^  -^0  Jle 

where  Me  is  the  mass  of  the  system  and  Ic  is  the  pitch  mass  moment  of  inertia  of  the  system 
about  the  mass  center  C. 

The  potential  energy  is  assumed  to  be  entirely  due  to  the  strain  energy.  In  view  of 
the  lay-up  symmetry,  it  is  reasonable  to  assume  that  uq  and  uq  are  zero  so  that  the  strain- 
displacement  relations  assume  the  simple  form 


7yz  =  V’y  +  ^,sf,  7x7  =  +  «l.x,  7xy  =  ^  {'^x,y  +  ^y,x)  e,  /) 

where  the  shear  strains  are  recognized  as  engineering  shear  strains.  Moreover,  the  symbols 
,  X  and  ,  y  in  the  subscripts  denote  derivatives  with  respect  to  x  and  y ,  respectively. 

For  the  jth  generally  orthotropic  layer  with  principal  axes  1  and  2,  the  constitutive 


equations  take  the  form 


■  0-1  ' 

J 

■  ei  ■ 

02 

C2 

T23 

=  [QY 

723 

Tl3 

713 

.'ri2. 

.712. 

(5) 


4 


The  elements  of  [Q]^’  are  related  to  material  properties  of  the  jth  layer  by 

M 


Oil  — 


0^2  — 


1  - 

Ei 


j  _  ^12-^2 

1  -  ^2*^1 

,  qL  =  g{, 


1  -  H;2*^i 
qL  =oL  Qie  =  ^i2 


(6a,  b) 

(6c,  d) 
(6e,/) 


where  Ej,  Gim  and  Uim  are  Young’s  moduli,  shear  moduli  and  Poisson’s  ratios,  respectively. 
The  stress,  engineering  strain  and  constitutive  relationships  for  the  ^th  layer  can  be  written 


as  follows: 


<t'  =  [ry]-v;',  c;-  =  [r.j-'e;',  aj-  =  [g.-ic;- 


(7a,  6,  c) 


where 


cos^  6j 

sin^  dj 

0 

0 

sin26y 

sin^  6j 

cos^  9j 

0 

0 

—  sin  29 j 

0 

0 

cos  9j 

-  sin  9j 

0 

0 

0 

sin  9j 

cos  9j 

0 

J  sin  29 j 

isin25j 

0 

0 

cos29j  . 

i!  =  K)-‘  [13,1  Kl 

-T 

ryz 

‘^xz 

^xy] 

T  _w 

7v* 

7ry 

7xy 

C"  : 

2 

2 

2  J 

’  '"j 

=  [o-l  (72  T23  Ti3 


ei  62 


723  713  712 

2  2  2  J 


(8a) 

(86) 
(8c,  (i) 
(8e,  /) 


in  which  the  conversion  to  engineering  strain  was  taken  into  account.  The  angle  9j  is  from 


the  plate  axis  x  to  the  material  axis  (l)y. 

Next,  we  propose  to  derive  the  laminate  streiin  energy.  The  laminate  consists  of  2k 
symmetrically-stacked  layers  of  variable  thickness.  Each  layer  height,  tj,  is  a  continuous 
function  of  i  and  y  (Fig.  2).  The  strain  energy  density  for  a  single  layer  j  or  its  counterpart 


below  the  midplane  is 


^  i=:l 

where  i  =  x,y,yz,  iz,  xy.  Summation  over  all  the  layers  and  integration  over  the  domain  of 

the  plate  leads  to  the  total  strain  energy 

1 


(10) 


5 


where  5  is  the  wing  span,  LE  =  LE{y)  the  wing  leading  edge  and  TE  =  TE{y)  the  wing 
trailing  edge.  Substituting  the  strain-displacement  relations,  Eqs.  (4),  and  the  constitutive 
relations,  Eq.  (7c),  into  Eq.  (10)  and  integrating  over  the  layer  thicknesses,  we  obtain  the 
total  strain  energy  expression 


+  (V’x.y  +  ^66  +  2t/>r,xV’y,y^2  +  (V»i,y  +  V’y.x)  ^6 

+  2t/)y,y  (t/)x,y  +  V’y.X )  ^5]  +  {tj  “  tj-l)  [(V’x  +  U^.x)^  ^5 

+  ii’y  +  ^44  +  2  (t/>,  +  to,*)  (V>y  + 


(11) 


It  must  be  recognized  here  that  tj  =  tj{x,y).  The  summation  can  be  eliminated  from  the 
strain  energy  expression  by  considering  the  total  laminate  extensional  and  bending  stiffness 
coefficients  and  Dab  defined  as 

=  t  {t,  -  C.1  =  5  E  (ij  -  t?-l)  ei.  (12) 

j=-k  ^  i=-fc 

and  we  observe  that,  because  the  thickness  of  the  various  layers  is  variable,  Aab  =  Aai(x,j/) 
and  Dab  =  Dab{x,  y).  Using  Eqs.  (12),  the  total  laminate  strain  energy  can  now  be  expressed 

in  the  compant  form 

V  =]r  I  I  (fijD^D  +  'Pa^^a)  dxdy  (13) 

JO  J LB 


where 


'Pd  = 


'Py,y 

,V’x,y  +  V’y.x. 


'Pa  = 


rpy+W^y' 

'Px  +  . 


(14a,  6) 


D  = 


1^11  Dn  Die 
Di2  D22  D26  1 
Die  D26  Dee. 


A45' 

>155. 


(14c,  d) 


It  remains  to  derive  an  expression  for  8W nc-  Under  consideration  is  a  wing  in  the  form 


of  a  trapezoidal  planform,  and  in  particular  one  characterized  by  low  aspect  ratio  and/or  for¬ 
ward  swept  configuration.  The  most  important  speed  regime  for  such  a  wing  is  undoubtedly 
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supersonic.  The  main  reason  for  including  the  aerodynamics  is  to  demonstrate  the  usefubess 
of  the  structural  model.  A  complete  mvestigation  of  a  wmg  would  require  appropriate  aero¬ 
dynamic  theories  for  subsonic,  transonic,  supersonic  and  perhaps  hypersonic  speed  regimes. 
The  usefubess  of  tbs  model  can  be  demonstrated  with  supersonic  aerodynamics,  chosen  for 
relevancy  and  for  relative  ease  of  application.  For  large  Mach  numbers,  »  1,  there 
is  a  weak  memory  effect,  in  addition  to  weak  tbee-dimensional  effects.  Tbs  opens  up  the 
prospect  of  a  point-function  relation  between  the  pressure  difference  p«  -  pi  and  the  dis¬ 
placement  w  of  the  wing,  wbch  makes  it  both  convement  and  usebl.  The  structural  model 
is  suitable  for  wings  of  any  aspect  ratio.  Wble  piston  theory  is  used  to  demonstrate  this 
model,  it  must  be  understood  that,  because  it  is  a  strip  theory,  it  is  not  as  suitable  for  very 
low  aspect  ratios  as  more  sophisticated  theories.  This  pomt  function  property  allows  the 
writing  of  an  explicit  a.erodynamic  operator  operating  on  the  vertical  displacement  of  the 
wing,  so  that  the  force  density  can  be  written  in  the  form® 


Ai,u.,  =  -c(x)  + 


where 


C(^)  =  ^ 


I  + 

2  dx 


(15) 


(16) 


in  which  M  is  the  Mach  number,  q  the  dynamic  pressure,  tfi  the  wing  half  thickness,  U  the 
free  stream  airspeed  and  7  the  ratio  of  specific  heats. 

The  nonconservative  virtual  work  is  due  to  the  aerodynamic  forces.  The  distributed  force, 
Eq.  (15),  multiplied  by  the  corresponding  virtual  displacement  yields  the  nonconservative 
virtual  work  density.  Then,  using  Eqs.  (Id)  and  (15)  and  integrating  over  the  domain  of  the 
wing,  we  obtab  the  nonconservative  virtual  work 

Te  /■* 


SWnc  =  f  I  AiiWpSwpdxdy  =  WcSwc  +  ^  WSwdxdy  (17) 

«/0  J LE 


where 


Wc  =  -  C'(x)|V’C  +  w,x  +  ^  +  {x-  xc)  '^dxdy  (18a) 
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is  a  resultant  aerodynamic  force, 


=  -  j*  C  (i)  |t/»c  +  UJ.X  +  ^  [^C  +  (a:  -  xc)  +  tb] }  (x  -  Xc)  dxdy  (185) 
is  a  resultant  aerodynamic  moment  about  an  axis  parallel  to  y  and  passing  through  C  and 


W  -  -C{x)  |v>c  +  to,*  +  ^  [the  +  ix-  xc)  ^c  +  ^]] 


(18c) 


is  an  aerodynamic  force  density. 

4.  THE  EIGENVALUE  PROBLEM 

Inserting  Eqs.  (3),  (13)  and  (16)  into  the  extended  Hamilton’s  principle  and  following  the 
usual  steps,®  we  can  obtain  the  boundary- value  problem  for  the  system.  The  boundary- value 
problem  consists  of  two  ordinary  differential  equations  for  u)o(t)  \j>c{t)  and  three  par¬ 
tial  differential  equations  for  to(i,y,t),  r/.*(i,y,t)  and  r/>j,(i,y,t),  together  with  appropriate 
boundary  conditions. 

The  differential  eigenvalue  problem  can  be  obtained  from  the  boundary- value  problem 
by  assuming  that  the  displacements  vary  exponentially  with  time.  The  state  of  the  art  does 
not  permit  a  closed-form  solution  of  the  differential  eigenvalue  problem,  so  that  one  must 
be  content  with  an  approximate  solution,  which  amounts  to  spatial  discretization  of  the 
problem  through  a  series  expansion.  As  a  result,  the  original  differential  eigenvalue  problem 
is  replaced  by  an  algebraic  eigenvalue  problem.  It  turns  out  that,  in  deriving  the  algebraic 
eigenvalue  problem,  it  is  more  efficient  to  use  directly  the  extended  Hamilton’s  principle,  Eq. 
(2).  To  this  end,  we  assume  a  solution  of  the  form 


wc{t)  =  ?l(i)>  V’c(<)  =  92(<) 

(19(1,6) 

n+2 

(19c) 

w{x,y,t)  =  <f>i{x,y)qi{t) 

2n+2 

(19,i) 

V>x(i,y,<)=  Y.  ^•(a:,y)gi(*) 

i=n-h3 

3n+2 

(19c) 

‘tf)y{x,y,t)  =  Y 

t=2ii+3 
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where  (f)i  axe  spa^e-dependent  trial  functions  and  g,-  axe  time-dependent  generalized  coordi¬ 
nates.  Inserting  Eqs.  (19)  into  Eq.  (3)  and  omitting  the  integral  limits  for  brevity,  we  can 
write  the  kinetic  energy  in  the  discretized  form 

7=1  (Afctbc  +  Me)  JJ mwdxdy  +  ^c  JJ ^c)  wdxdy 

+  ^II  rnw^dxdy  =  ^q^Mq  (20) 

where  q  is  a  3n-vector  of  amplitudes  and  M  is  the  synunetric  mass  matrix  having  the 
components 


Mil  =^Ci  ^12  =  0 

(21.1,4) 

Mij  =  jj  TTKfijdxdy,  3  <  i  <  n  +  2 

(21c) 

Miy  =.  0,  i  >  n  -t-  2 

(21d) 

M21  =  0,  M22  =  Ic 

(21e,/) 

M2j  =  JJ  m(i  -  xc)<f>}dxdy,  3  <  j  <  n  -f  2 

(21s) 

^2]  =  0)  j  >  n  +  2 

(214) 

Mij  =  JJ  m4)i4>jdxdy,Z  <i,j  <n  +  2 

(21i) 

Mij  =  0,  ij  >  n  +  2 

(21j) 

Moreover,  introducing  Eqs.  (19c-e)  into  Eq.  (13)  and  considering  Eqs. 

(14),  we  obtain  the 

discretized  potential  energy 

V=^JJ  (t/>St/)zj  +  ^aMa)  dxdy  = 

(22) 

where  K  is  the  symmetric  stiffness  matrix  having  the  entries 

Kij  =0,  i  =  1, 2  and  1  <  ;  <  3n  +  2 

(23a) 

Kij  =  JJ  [<f>i,x<f>j,xMs  +  {<f>i,y(i>j,x  +  <i>i,x<i>j,y)  Mh 

dxdy,  3  <  t,  j  <  n  -|-  2 

(236) 

Kij  =JJ  {<f>i,x4>jMs  +  <i>MjMh)dxdy, 
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(23c) 


3<t<n  +  2,  n  +  3<j<2n  +  2 

Kij  =  JJ 

3<t<n  +  2,  2n  +  3<  j  <3n  +  2  (23(i) 

Kij  =  JJ  [<l>i<j>jAss  +  (j)i,x(f>j,xDli  +  <j>i,y<f>j,yDs6 

+  +  4i,x4>i,y)  -^>16]  n  +  3  <  t,;  <  2n  +  2  (23e) 

Kij  =  JJ  {<f)i4)jMi  +  (i>i,x<f>j,xDu  +  4>i,y<i>j,xDs6  +  4>i,x<i> j,y Dl2 

+4>i^y<f>j,yD2t)  dxdy,  n  +  3  <  i  <  2n  +  2,  2n  +  3  <  J  <  3n  +  2  (23/) 

Kij  =  JJ  [<f>i<l>jA^4  +  4>i,x<l>j,xD66  +  <i>i,y<i>j,yD22 

+  +  4>i,x<f>j,y)  -^26]  dxdy,  2n  +  3  <  ij  <  3n  +  2  (23y) 

Finally,  inserting  Eqs.  (19c-  e)  into  Eq.  (17)  and  considering  Eqs.  (18),  we  can  write 
the  discretized  virtual  work  due  to  aerodynamic  forces  in  the  form 

Wr,c  =  WcSwc  +  +  JJ  WSwdxdy  =  -6q^  (^Ka^  +  (24) 

where  the  entries  of  the  nonsymmetric  matrix  Ka  are 

KAij  =  0,  ;  =  l,2,...,n  +  2  (25a) 

Ka21  =  JJ  Cdxdy  (256) 

ka22  =  JJ 

KA2j  =  JJ  C<f>jdxdy,  3  <  j  <  n  +  2  (25d) 

Kah  =  // C<l>i,xdxdy,  d<i<n  +  2  (25e) 

KAi2  =  JJ  <f>i,xdxdy,  3<i<n  +  2  (25/) 

KAij  =  JJ  C<f>i,r4>jdxdy,  2<iJ<n  +  2  (25y) 

KaH  =  0,  1  <  i  <  n  +  2,  j  >  n  +  2  (25/i) 

KAij  =  0,  i>n  +  2  (25i) 

and  those  of  the  symmetric  matrix  H  are 

Hii  =  JJcdxdy,  Hn  =  JJc  {x  -  xc)  dxdy  (26a,  6) 
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H\j  =  Jj  C<f)jdxdy,  3  <  j  <  n  +  2 

(26c) 

H22  =  jJ  C{x-xcf  d.xdy 

(26<i) 

H2j  =  //^  (^  ”  ®^)  3  <  j  <  n  +  2 

(26e) 

Hij  =  JJ  C(f)i(f>jdxdy,  3  <  t,;  <  n  -f  2 

(26/) 

Hij=0,  1  <  i  <  n  +  2,  j  >  n  +  2 

(265) 

Hij  =0,  t  >  n  -(-  2 

(26/i) 

Introducing  Eqs.  (20),  (22)  and  (24)  into  Eq.  (2),  integrating  by  parts  with  respect  to 
time  and  invoking  the  arbitrariness  of  (Jq,  we  obtain  a  set  of  simultaneous  ordinary  differential 

equations,  which  can  be  written  in  the  compact  form 

Mq  4-  +  {K  +  Ka)  q  =  0 

(27) 

Then,  letting 

q{t)  =  ae^* 

(28) 

we  obtain  the  desired  algebraic  eigenvalue  problem 

(x^M  -i-  —H  +  K  +  Ka^  a  =  0 

(29) 

and  we  note  that  the  problem  is  nonself-adjoint  due  to  the  presence  of  the  aerodynamic 
matrices  H  and  Ka- 

The  solution  of  the  eigenvalue  problem  (29)  requires  that  it  be  cast  in  state  form.  To  this 
end,  we  introduce  the  state  vector  x  =  [a^  Aa^]^.  Then,  adjoining  the  identity  Aa  =  Aa, 

Eq.  (29)  can  be  rewritten  in  the  state  form 

Ax  =  XBx  (30) 


where 


A  = 


0 

-{K  +  Ka) 


I  ■ 
-H/U. 


B  = 


I  0 
0  M 


(31o,6) 
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Because  the  matrix  B  is  singular,  some  of  the  eigenvalues  are  infinite.  Numerical  solutions 
for  this  type  of  eigenvalue  problems  can  be  obtained  by  an  algorithm  described  in  Ref.  10. 

5  approximate  solutions  in  terms  of  QUASI-COMPARISON 
FUNCTIONS 

In  Sec.  4,  we  derived  a  discretized  model  for  our  system  by  assuming  an  approximate 
solution  in  the  form  of  a  linear  combination  of  trial  functions.  The  accuracy  of  the  approx¬ 
imate  solution  and  the  computational  time  depend  on  the  nature  of  the  trial  functions.  To 
examine  the  nature  of  the  trial  functions,  a  brief  discussion  of  the  various  classes  of  functions 
should  prove  beneficial. 

Exact,  closed-form  solutions  of  differential  eigenvalue  problems  represent  the  class  of 
eigenfunctions.  Clearly,  they  satisfy  both  the  differential  equations  and  all  the  boundary 
conditions  of  the  problem.  In  most  practical  problems  closed-form  solutions  do  not  exist. 
They  certainly  do  not  exist  in  the  problem  at  hand.  Approximate  solutions  generally  fail 
to  satisfy  the  differential  equations.  Functions  satisfying  certain  differentiability  conditions 
and  all  the  boundary  conditions,  but  not  necessarily  the  differential  equations,  are  called 
comparison  functions}  In  many  problems,  including  that  under  consideration,  comparison 
functions  are  diflBcult  to  generate.  In  the  case  of  self-adjoint  systems,  it  is  often  advantageous 
to  formulate  the  eigenvalue  problem  by  a  vaxiational  approach,  which  amounts  to  rendering 
Rayleigh’s  quotient  stationary.  In  this  case,  the  trial  functions  need  satisfy  reduced  differen¬ 
tiability  conditions  and  the  geometric  boundary  conditions  alone.  Such  functions  are  known 
as  admissible  functions}  It  has  been  demonstrated  in  Ref.  7  that  solutions  in  terms  of  ad¬ 
missible  functions  can  at  times  converge  very  slowly.  To  improve  convergence,  a  new  class  of 
functions  was  created  in  Ref.  11,  namely  the  class  of  quasi-comparison  functions.  The  quasi¬ 
comparison  functions  are  linear  combinations  of  admissible  functions  capable  of  satisfying 
all  the  boundary  conditions  of  the  problem.  Quasi-comparison  functions  can  be  generated 
by  combining  several  families  of  admissible  functions  chosen  so  as  to  permit  the  satisfaction 
of  the  natural  boundary  conditions.  It  turns  out  that  solutions  in  terms  of  quasi-comparison 
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functions  tend  to  satisfy  not  only  the  natural  boundary  conditions  but  also  the  differential 
equations  much  more  accurately  than  solutions  in  terms  of  mere  aximissible  functions.  In 
fact,  solutions  in  terms  of  quasi-comparison  functions  can  at  times  be  more  accurate  than 
solutions  in  terms  of  comparison  functions,  because  combinations  of  functions  from  several 
families  with  different  characteristics  permit  better  approximations  of  the  solution  through¬ 
out  the  interior  of  the  domain  than  combinations  of  functions  from  a  single  family.^-“  In 
Refs.  12  and  13,  it  was  demonstrated  that  solutions  in  terms  of  quasi-comparison  functions 
exhibit  the  same  superior  convergence  characteristics  in  the  case  of  nonself-adjoint  systems 

as  well. 

In  view  of  the  above  discussion,  we  propose  to  choose  the  trial  functions  in  Eqs.  (19c-  e) 
from  the  class  of  quasi-comparison  functions.  These  functions  must  depend  on  both  i  and  y 
and  are  assumed  to  have  the  form  of  products  of  chordwise  functions  Xkj{x)  and  spanwise 

functions  Yi^{y). 

The  eigenvalue  problem  is  defined  in  terms  of  trial  functions  <^i  and  <i>j  (i,  j  =  1,2,3, . 
3;^  +  2),  where  n  is  the  number  of  trial  functions  for  each  of  the  three  displacements  ui,  V’l 
and  i/jy.  There  are  additionally  two  rigid-body  modes,  namely,  <j>i  =  I  and  <^2  =  1- 

Let  the  same  set  of  trial  functions  be  used  for  each  of  the  displacements  u;,  V**  and  0,,. 
The  undetermined  coefficients,  of  course,  are  not  the  same.  Then,  for  i  =  3, 4, 5, . . . ,  3n  -H  2, 

<f>f{x,y)  =  Xkf  {x)Yej  (y) ,  l<  f  <n 

where 

/  =  i-2-n.INT(i^) 

Similarly,  for  the  companion  set  of  trial  functions 

=  Xk,{x)Yt^{y).  1  <  5  <  n 

in  which 

for  g  =  j  —  2-n-  INT  ^  ^ 


(326) 

(32c) 

(32d) 
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where  INT(  )  implies  truncation  to  the  corresponding  integer  value. 

The  functions  Xkj  and  Yi^  are  each  chosen  from  two  different  families.  The  indices  kf 
and  if  determine  the  order  of  combination  of  functions  and  into  and  hence  they 
depend  on  /,  which  in  turn  depends  on  i  according  to  Eq.  (32b).  Because  the  interest  is 
in  low-aspect  ratio  wings,  the  dependence  of  kf  and  if  on  /  is  such  that  the  spanwise  and 


chordwise  functions  are  about  equal  in  number  or 

=  /  =  1,2,3,...,^  (33a) 

if  =  2  —  n  —  (336) 

kf  =  l-d-if  /=l,2,3,...,n  (33c) 

where 

d  =  NINTy^  /=  1,2,3,...,^  (34a) 

d  =  NINTy2/  -  n  /=  ^  +  1,  ^  +  2,....,n  (346) 


and  NINT(  )  implies  rounding  to  the  nearest  integer  value.  It  is  assumed  in  Eqs.  (33)  and 
(34)  that  n  is  an  even  integer. 

The  chordwise  functions  Xkj{x)  are  segments  of  a  sine  function  with  an  appropriate 
number  of  zero  crossings.  Accounting  for  wing  sweep  and  the  variable  limits  of  i,  X^j 
actually  become  functions  of  both  i  and  y.  The  first  family  of  chordwise  function  are  chosen 


cLS 


Xkj{x,y)  = 


X  -  y  tan  rn 


I - ^ LE<x<TE  (35) 

[r  -  y(tan77i  -  tan  tit)}  J 

The  spanwise  functions  Yt^{y)  are  also  appropriate  segments  of  a  sine  function  with  an 

appropriate  number  of  zero  crossings,  or 


=  -  sin 


{2if  -  l)7r 


23 


,0<y<3 


(36) 


The  second  family  of  functions  Xkj  consists  of  terms  from  a  power  series,  except  that  the 
functions  alternate  direction  between  the  leading  edge  and  the  trailing  edge.  The  wing  sweep 
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requires  shifting  and  scaling  as  before,  resulting  in 

Xi/(x,y)  =({r  [l  +  (-1)*/]  -  (-1)*^^  +  [l  +  (-1)*^]  -  ^^VL)y}/2[r 

+  (tanr/r  -  tanr;^)  y]  )  ^  \  LE<x<TE  (37) 


where  INT(  )  implies  truncation  to  the  corresponding  integer  value.  Equation  (37)  appears 
complicated  but  it  merely  represents  a  power  series  shifted,  scaled,  and  alternating  in  direc¬ 
tion  so  as  to  accomodate  the  swept  boundaries.  The  second  family  of  spanwise  functions  is 

simply 

0<!,<3  (38) 

The  algebraic  eigenvalue  problem  is  assembled  inside  a  double  loop  for  /  and  g  from  1 
to  n.  The  entries  of  the  matrices  M,  K,  Ka  and  H  are  given  by  Eqs.  (21,)  (23),  (25)  and 


(26),  respectively. 

The  complete  set  of  two-dimensional  quasi-comparison  functions  is  expressed  in  terms 
of  /  through  if  and  kf  as  follows; 


V  ^  "b  2  .  - 

<f>f  sin  ay,  /  =  3,4,...,  «/ =  1 

(f)f  sin  ay,  /  =  3,4,...,—^,  kf  =  2 
^  rs  ^ 

y.  .b  +  cx  +  dy  n  +  2  r 

(if  =- sin  ay  sin - ,  /-3,4,...,  ,  «/ £  o 

^  s  r  +  ey  ^ 


(39a) 


(39c) 


where 


(40a) 

(406) 

. ,  n  ,  kf  >Z 

(40c) 

,  3r7r 

b  = - ,  e  =  tan  tjt  -  tan  r]L 

4 

(41a,  6) 

(2if  -  l)  7r 

p_-tanJ7i,  a- 

(41c,  d) 
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(41c,/) 


c^^(kf-l),d  =  ^[le+{k,-l)p 

(3  =  -  (-Lf> ,  s  =  INT 

Similar  series  can  be  written  for  (f>g  by  replacing  /  by  and  i  by  j.  The  simpler  forms 
for  <f>f  when  A:/  =  1,2  are  due  to  the  simple  form  of  the  rigid  chordwise  shape  functions 
corresponding  to  the  first  two  terms  in  Xkf{x),  as  mentioned  previously.  An  example  of  the 
series  (f>f{x,  y)  for  n  =  6  should  prove  helpful.  The  six  ternas  consist  of  three  terms  from  each 
of  the  two  families,  or 


(42a,  6,  c) 
(42ti,  e,  /) 


6.  NUMERICAL  RESULTS 

In  presenting  numerical  results,  it  is  useful  to  define  several  dimensionless  quantities. 
The  frequency  can  be  nondimensionalized  as  follows: 

Similaxly,  a  speed  paxameter  can  be  defined,  yielding  the  dimensionless  dynamic  pressure 

2qA^  ,44) 

The  quantities  (  )ref  reference  values.  The  quantities  and  pjef  axe  E2  and  p  for 
the  main  structural  material,  t^ef  is  the  total  thickness  at  the  wing  root  and  leading  edge 
intersection.  Other  dimensionless  parameters  axe  taper  ratio  TR  and  the  aspect  ratio  AR. 
Extensive  vibration  and  flutter  results  can  and  have  been  generated  with  this  model.  The 
intent  here  is  to  develop  and  present  the  formulation,  so  that  the  presented  results  axe  merely 


a  summary. 
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i.  Convergence  and  accuracy 

The  convergence  of  the  free- vibration  eigenvalues  of  an  anisotropic  low  aspect  ratio  wing 
is  considered  in  Fig.  3.  The  convergence  rate  of  the  first  five  nonzero  eigenvalues  is  depicted 
as  the  frequency  parameter  fl  versus  the  number  of  terms  in  the  approximating  series.  The 
first  two  eigenvalues  have  zero  value  and  correspond  to  rigid-body  modes,  which  converge 
immediately;  they  are  not  shown  in  Fig.  3.  It  should  be  pointed  out  here  that  a  Hneax 
combination  of  admissible  functions  must  possess  a  minimum  number  of  terms  before  it 
qualifies  as  a  quasi-comparison  function.  The  number  must  be  such  as  to  permit  satisfaction 
of  all  natural  boundary  conditions.  As  soon  as  this  number  of  terms  has  been  reached, 
convergence  of  the  computed  eigenvalues  to  the  actual  ones  is  relatively  rapid. 

The  two  flutter  mechanisms  of  most  concern  are  the  bending-torsion  mode  and  the  body- 
freedom  mode.  The  body-freedom  flutter  mode  occurs  instead  of  the  divergence  mode  when 
the  pitch  degree  of  freedom  is  introduced.  Reasonably  accurate  representations  of  the  pitch 
mode  and  the  first  two  flexible  modes,  as  well  as  fairly  accurate  representation  of  the  third 
flexible  mode,  are  required  to  capture  these  two  mechanisms.  The  indications  from  Fig.  3 
are  that  20  terms  are  sufficient  for  this  purpose,  certainly  for  establishing  a  trend. 

Convergence  of  the  flutter  dynamic  pressure  parameter  for  a  forward  swept  anisotropic 
wing  is  shown  in  Fig.  4.  The  flutter  mechanism  acting  is  body-freedom  flutter.  Reasonable 
convergence  has  been  achieved  with  20  terms. 

Figure  5  shows  three  frequencies  corresponding  to  an  isotropic  plate  of  aspect  ratio  5; 
the  computed  values  are  compared  with  the  exact  values  given  in  Table  11.4  of  Ref.  14. 
The  plate  is  1  /2  inch  steel  and  is  fixed  on  a  short  side.  The  remaining  three  edges  are  free. 
The  clamped  condition  is  simulated  with  the  present  model  by  taking  very  large  values  for 
fuselage  mass  and  pitch  inertia.  The  results  are  quite  acceptable.  As  expected  from  the 
discussion  of  quasi- comparison  functions  the  matching  is  poor  for  small  numbers  of  terms 
but  begins  to  converge  quite  rapidly  when  a  certain  critical  number  of  terms  is  reached. 

The  accuracy  of  the  flutter  analysis  is  checked  through  a  comparison  with  results  obtained 
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in  Ref.  15  by  Rossettos  and  Tong,  which  is  one  of  the  few  works  available  on  the  flutter 
of  a  plate  with  free  edges  in  the  supersonic  region.  Their  model  consists  of  a  single  layer, 
constant  thickness,  square  composite  plate  clamped  at  one  edge  and  subjected  to  supersonic 
flow  on  the  upper  surface.  The  aerodynamics  is  based  on  a  variant  of  the  piston  theory  that 
converges  to  Eq.  (15)  for  sufficiently  high  Mach  numbers  and  with  flow  over  both  surfaces. 
Comparison  is  for  a  square  composite  plate  with  a  fiber  orientation  of  24°.  The  model  derived 
here  predicted  on  the  average  flutter  speeds  6%  higher  than  those  of  Ref.  15. 

ii.  Free  vibration  results 

The  interest  in  the  free  vibration  ca^e  is  mainly  for  its  relation  to  the  flutter  results. 
The  two  flutter  mechanisms,  body-freedom  and  bending-torsion,  dictate  that  concentration 
is  on  the  pitch  mode  and  the  first  two  flexible  modes,  although  the  first  few  flexible  modes 
above  these  will  also  have  some  influence.  The  convergence  of  the  rigid-body  modes,  pitch 
and  plunge,  is  not  an  issue. 

For  body-freedom  flutter,  the  most  critical  parameter  is  geometrical.  If  the  wing  is  not 
swept  forward  or  nearly  so  then  this  flutter  mechanism  cannot  arise.  Work  by  a  number 
of  authors  demonstrate  a  surprising  influence  of  a  tailoring  layer  swept  just  a  few  degrees 
off  nominal.  Body-freedom  flutter  can  be  eliminated,  or  pushed  up  to  significantly  higher 
speeds,  if  a  percentage  of  the  spanwise  plies  are  reoriented  10°  forward.  Two  factors  are  at 
work  in  this  case  and  both  are  free  vibration  matters.  When  an  orthotropic  layer  for  which 
E\»  E2  is  reoriented,  bending  begins  to  occur  in  the  direction  normal  to  Ei.  The  result 
is  an  effective  rotation  of  the  elastic  axis  in  the  direction  opposite  to  the  reorientation  of 
the  tailoring  ply.  The  other  factor  at  work  is  the  necessary  reduction  of  spanwise  bending 
stiffness  that  must  accompany  the  reorienting  of  plies  away  from  the  midchord.  This  second 
factor  eventually  overpowers  the  first,  so  that  ply  reorientations  are  most  effective  when  they 
are  small. 

Figure  6  shows  the  frequency  parameter  of  the  first  flexible  mode  varying  with  the  sweep 
angle  of  a  tailoring  layer.  The  results  for  three  wings  are  shown.  An  interesting  phenomenon 
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is  observed  in  Fig.  6  for  the  low-aspect  ratio  wing.  This  wing  has  its  highest  bending  stiffness 
associated  with  a  significantly  aft  swept  tailoring  layer,  «  =  10”.  The  implication  is  that 
tailoring  phes  might  not  be  as  effective  for  low-aspect  ratio  wings.  This  trend,  shown  to  be 
true  in  the  next  subsection,  was  not  previously  discovered  due  to  inherent  limitations  of  the 
one-dimensional  models  commonly  used. 

Much  has  been  said  about  the  benefits  of  using  the  anisotropicity  of  a  composite  layup 
to  tailor  the  response  of  a  wing.  Such  discussion  has  been  limited  almost  entirely  to  ply  ori¬ 
entation.  The  dominant  use  of  one-dimensional  structural  models  nas  forced  this  limitation 
on  the  current  investigation.  One  can  easily  suppose  that  distribution  of  the  tailoring  ma¬ 
terial  in  both  the  spanwise  and  chordwise  directions  might  significantly  affect  the  response. 
These  trends  are  readily  investigated  with  the  current  plate  model.  As  ply  distribution  for  a 
low-aspect  ratio  anisotropic  wing  is  varied  from  a  concentration  near  the  traiUng  edge  to  the 
leading  edge,  there  is  a  change  in  frequency  of  the  natural  modes.  The  associated  frequency 
paxameters  increase  by  8%,  6%  and  3%,  respectively,  for  the  first  three  flexible  modes.  The. 
bending-torsion  separation  increases  by  5%.  For  a  similar  spanwise  variation  from  wing  tip 
concentration  to  wing  root,  frequencies  again  vary.  In  this  case,  the  increases  are  34%,  14% 
and  12%,  respectively,  for  the  first  three  flexible  modes.  For  bending-torsion  separation,  the 
change  is  slight  at  2.5%.  This  reflects  the  increase  in  bending  stiffness,  which  occurs  when 
material  is  shifted  near  the  wing  root.  Clearly,  the  two-dimensional  nature  of  this  material 
distribution  has  an  effect  on  the  free  vibration  outcome.  It  seems  reasonable  to  expect  that 
the  flutter  response  will  also  be  affected,  aUowing  a  refinement  of  the  concept  of  tailoring 

with  composites. 

iii.  Flutter  analysis  results 

Two  flutter  mechanisms  axe  of  primary  concern  for  a  symmetric  model  of  a  wing  with 
forward  or  aft  sweep.  The  first  is  the  classical  bending- torsion  flutter  that  can  occur  for  a 
wing  of  any  sweep  when  the  frequencies  of  these  two  modes  begin  to  coalesce  as  the  speed 
increases.  One  mode,  the  torsion  mode,  turns  downward  toward  a  very  stable  condition  just 
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as  the  bending  mode  heads  for  instability.  The  two  frequencies  coalesce  as  the  instability 
approaches.  Complete  coalescence  does  not  occur  when  the  fuselage  is  capable  of  rigid-body 
motions^  The  rigid-body  pitch  frequency  assumes  a  nonzero  value,  which  rises  as  the  speed 
increases.  Without  a  tailored  layup,  the  pitch  and  bending  modes  would  surely  combme 
to  cause  body-freedom  flutter  before  the  bending-torsion  coupling  could  arise.  Recall  that 
bending-torsion  flutter  is  unlikely  at  any  speed  for  a  significant  forward  sweep  angle.  It 
is  seen  then  that  the  same  tailoring  which  eliminates  body-freedom  flutter  by  apparent  aft 
sweeping  of  the  elastic  axis  also  makes  bending-torsion  flutter  possible.  A  trade-off  is  implied 
and  caution  must  be  exercised  when  contemplating  a  tailored  composite  design  solution. 

Next,  the  characteristics  of  body-freedom  flutter  are  considered.  For  a  swept  forward 
orthotropic  wing,  the  flutter  mode  wiU  be  body-freedom  flutter,  or  divergence  if  clamped. 
The  bending-torsion  case  occurs  at  a  much  higher  speed  or  may  not  occur  at  all  for  a  large 
forward  sweep  angle.  The  pitch  mode  is  characteristically  very  lightly  damped  and  will 
move  toward  instability  as  the  bending  mode  moves  away  from  the  imaginary  axis.  The 
pitch  and  bending  frequencies  coalesce  eventually  with  increasing  dynamic  pressure  rather 
than  separating,  as  would  be  the  case  for  bending-torsion  flutter. 

Weisshaar  et  al.^  have  reported  a  trend  in  the  ratio  fiflutter/I^fundamental  for  the  two  flutter 
mechanisms  just  presented.  In  a  particular  case  in  Rnf.  5,  the  ratio  is  equal  to  0.26  for  a 
body-freedom  flutter  and  to  3.40  for  a  bending-torsion  flutter.  The  ratios  for  comparable 
cases  in  this  work  are  0.42  and  2.64,  respectively. 

The  effect  of  tailoring  plies,  particularly  on  divergence  or  body-freedom  flutter  has  re¬ 
ceived  a  great  deal  of  attention.  Figure  7  shows  the  flutter  dynamic  pressure  versus  the 
tailoring  ply  angle  for  two  orthotropic  wings  {AR  =  3  and  AR  =  6)  and  one  isotropic  wing 
{AR  =  8).  In  all  three  cases,  body- freedom  flutter  is  critical  and  small  negative  ply  angles 
are  expected  to  be  effective.  There  is  an  implication  in  Fig.  6  that  the  low-aspect  ratio 
wing  does  not  realize  much  benefit,  because  its  maximum  bending  stiffness  occurs  for  a 
somewhat  more  aft  swept  configuration  than  the  other  higher-aspect  ratio  cases.  This  fact. 
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demonstrated  in  Fig.  7,  has  important  design  impUcations  for  the  design  of  low-aspect  ratio 
wings.  If  the  model  used  does  not  account  for  low-aspect  ratios,  then  this  phenomenon  wiU 

be  entirely  missed. 

Figure  7  shows  the  interesting  and  somewhat  surprising  result  that  flutter  dynamic 
pressure  is  strongly  affected  by  very  small  changes  in  the  layer  orientation  in  the  vicinity  of 
0°.  Virtually  all  the  analytical  and  experimental  work  on  tailoring  forward  swept  wings  to 

date  does  not  consider  low-aspect  ratio  wings. 

Finally,  the  distribution  of  tailoring  plies  over  the  wing  planform  is  considered.  The 
present  model  is  ideaUy  suited  to  investigation  of  parameters  such  as  layer  thicknesses,  which 
are  free  to  vary  in  both  chordwise  and  spanwise  directions.  Such  distribution  investigations 
are  likely  to  result  in  lower  structural  weight  or  higher  flutter  speeds,  so  that  they  must  be 
included  in  any  search  for  an  optimal  wing  design.  The  present  model  can  be  used  to  conduct 
such  investigations  with  ease.  The  results  are  impressive.  In  the  chordwise  direction,  a  flutter 
dynamic  pressure  increase  of  25%  or  decrease  of  25%  is  observed,  depending  on  whether  the 
leading  or  traiUng  edge  is  favored  for  a  linear  material  distribution,  respectively.  Similarly, 
in  the  spanwise  direction,  the  flutter  dynamic  pressure  increases  29%  if  the  wing  root  is 
favored,  but  drops  35%  if  the  plies  are  concentrated  near  the  tip.  These  results  are  for  an 
orthotropic  wing  of  aspect  ratio  6  and  a  forward  sweep  angle  of  20°.  The  tailoring  layer  is 
swept  10°  forward  of  the  midchord  line.  These  new  results  go  a  long  way  toward  justifying 

this  more  sophisticated  model. 

7.  SUMMARY  AND  CONCLUSIONS 

At  the  outset,  the  goal  was  to  develop  a  structural  model  representing  a  modern  aircraft 
wing  sufficiently  well  and  one  that  is  sufficiently  simple  to  be  used  in  a  multidiscipHnary 
optimization  problem.  This  goal  has  been  accompUshed  with  the  development  of  a  plate 
model  that  accounts  for  shear  deformation,  variable  layer  thickness  and  orientation,  fuselage 
degrees  of  freedom  and  chordwise  flexibility.  Consistent  with  conclusions  reached  by  previous 
investigators,  certain  factors  such  as  wing  camber,  unsymmetric  ply  stacking  and  rotatory 
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inertia  are  sufficiently  small  for  the  problem  at  hand  that  they  can  be  ignored. 

The  algebraic  eigenvalue  problem  was  formulated  directly  from  the  extended  Hamil¬ 
ton  principle.  In  solving  the  algebraic  eigenvalue  problem,  the  recently  developed  quasi¬ 
comparison  functions  were  relied  upon  to  guarantee  rapid  convergence. 

The  numerical  results  obtained  by  means  of  the  current  model  confirmed  various  trends 
established  by  earlier  investigators.  They  demonstrate  that  there  is  a  very  wide  array  of  inter¬ 
related  parameters  affecting  the  instability  speed.  The  most  clear  outcome  is  a  demonstration 
of  the  necessity  of  a  formal  optimization  approach  to  bring  consistency  to  the  investigation 
of  all  the  pertinent  parameters. 

The  numerical  results  are  limited  to  demonstration  of  the  various  features  of  the  model 
and  opening  forays  into  regions  of  investigation  previously  closed  to  simpler  models.  Nev¬ 
ertheless,  several  previously  unreported  trends  having  significant  bearing  on  the  design  of 
modern  wings  were  revealed.  These  initial  inquiries  open  the  door  to  further  research  and 
refinement  of  the  concept  of  composite  tailoring.  The  results  indicate  that  composite  tailor¬ 
ing  of  forward  swept  wings  may  not  be  nearly  as  effective  for  low-aspect  ratio  wings  as  for 
high-aspect  ratio  wings.  This  is  an  important  result  which  virtually  mandates  a  structural 
model  of  at  least  the  level  of  sophistication  used  here  to  accomplish  optimal  design  of  a 
low-aspect  ratio  forward  swept  wing. 

A  second  and  also  important  new  finding  pertaining  to  two-dimensional  models  is  the 
merit  of  tailoring  ply  distributions.  Beyond  finding  the  optimum  angle  of  orientation  of 
a  tailoring  ply,  significant  improvements  in  the  flutter  speed  are  left  undiscovered  if  the 
distribution  is  not  also  subjected  to  scrutiny.  The  results  presented  here  indicate  that  the 
tailoring  plies  are  most  effective  when  the  distribution  favors  the  leading  edge  and  the  wing 
root.  This  result  is  also  unobtainable  with  more  conventional  models  and  points  once  again 
to  the  need  of  sufficient  sophistication  in  modeling.  Optimization  methods  are  dictated  by 
a  large  number  of  interacting  parameters,  but  such  effort  is  wasted  if  the  physical  system  is 
not  modeled  with  sufficient  accuracy. 
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The  integrated  design  of  a  structure  and  its  control  system  is  treated  as  a  multiobjective  optimization  problem. 
Structural  mass  and  a  quadratic  performance  index  constitute  the  vector  objective  function.  The  closed-loop  per¬ 
formance  index  is  taken  as  the  time  integral  of  the  Hamiltonian.  Constraints  on  natural  frequencies,  closed-loop 
damping,  and  actuator  forces  are  also  considered.  Derivatives  of  the  objective  and  constraint  functions  with 
respect  to  structural  and  control  design  variables  are  derived  for  a  finite  element  beam  model  of  the  structure 
and  constant  feedback  gains  determined  by  independent  modal  space  control.  Pareto  optimal  designs  generated 
for  a  simple  beam  demonstrate  the  benefit  of  solving  the  integrated  structural  and  control  optimization  problem. 


Introduction 

HE  fields  of  structural  optimization  and  optimal  control  have 
made  significant  strides  with  the  advent  of  the  computer  age. 
Each  discipline  has  independently  generated  useful  methods  and 
computational  tools  for  design.  Yet.  they  share  an  approach  in 
which  certain  system  parameters  or  design  variables  are  selected 
so  as  to  optimize  an  objective  function  subject  to  a  given  set  of 
constraints.  In  structural  optimization,  the  objective  function  is 
often  related  to  the  cost  of  the  structure,  and  it  may  involve  the 
\\  eight  of  the  structure  or  the  \oIume  of  material.  The  design  \'ari- 
ables  typically  characterize  the  material  distribution  or  geometry 
of  the  structure.  Any  quantity  characterizing  the  response  of  the 
structure,  such  as  stress,  displacement,  or  frequency,  may  be  con¬ 
strained  to  preclude  structural  failure.  The  structural  response  may 
be  static  or  dynamic,  although  time  usually  plays  no  panicular  role 
in  the  formulation  of  the  structural  optimization  problem.’-  In 
contrast,  in  modem  optimal  control  a  control  law  optimizing  some 
performance  index  (the  objective  function)  over  a  given  time  inter¬ 
val  is  synthesized.  In  the  linear  quadratic  regulator  (LQR)  theory,^ 
control  gains  relating  actuator  forces  to  sensor  outputs  by  means  of 
a  linear  transformation  are  typical  control  design  variables.  The 
performance  index  is  rendered  independent  of  time  by  integrating 
over  the  control  time. 

In  common  design  practice,  structural  design  precedes  control 
design.  The  mathematical  model  for  the  structural  optimum,  or  at 
least  the  final  structural  design,  constitutes  the  plant  for  the  control 
design.  The  control  designer  then  proceeds  to  synthesize  the  opti¬ 
mal  control  system  for  the  given  plant.  Uncertainties  in  structural 
modeling  of  the  plant  are  often  considered  to  ensure  that  the  con¬ 
trol  system  is  robust  with  respect  to  plant  errors.  The  possibility 
exists  that  the  plant  can  be  altered  so  as  to  enhance  performance  in 
conjunction  with  both  structural  and  control  optimization.  As  a 
result  of  difficulties  in  lifting  and  deploying  heavy  objects  such  as 
space  stations,  which  can  contain  large  solar  arrays,  antennas  and 
precision  laser,  or  optical  systems,  the  spacecraft  structures  must 
be  highly  flexible.  Moreover,  stringent  performance  requirements 
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for  pointing  accuracy,  vibration  suppression,  shape  control,  etc., 
demand  active  controls  to  augment  any  passive  damping.  The  goal 
of  integrated  design  is  to  take  advantage  of  any  synergistic  interac¬ 
tion  between  the  flexible  structure  and  its  active  control  system. 

The  integrated  structural  and  control  design  problem  was  first 
considered  "by  Hale  et  al."'  Messac  and  Turner  emphasized  trans¬ 
formation  of  the  state  space  to  modal  coordinates  and  the  depen¬ 
dence  of  the  weighting  matrices  on  the  modes.  In  formulating  the 
combined  optimization  problem  Salama  et  al.*  established  a  prece¬ 
dent  followed  by  many  others.  They  eliminated  control  design 
variables  by  considering  steady-state  (constant)  gains  and  select¬ 
ing  particular  weighting  matrices  (identity)  for  the  quadratic  per¬ 
formance  index. 

Haftka  et  al.’  minimized  control  effort  through  structural 
changes  while  maintaining  specified  damping  ratios.  However,  the 
damping  ratios  of  the  lower  modes  decreased,  and  the  frequency 
response  magnitude  and  decay  time  increased.  In  another  early 
study,  Venkayya  and  Tischler*  reached  a  different  conclusion, 
namely,  that  a  nominal  truss  model  and  a  structurally  optimal  truss 
had  nearly  the  same  control  design  for  regulating  a  disturbance. 
They  considered  a  disturbance  dominated  by  the  first  mode,  but 
both  structural  designs  had  the  same  fundamental  eigenvalue,  per¬ 
haps  explaining  the  similar  control  designs.  Miller  and  Shim?  also 
considered  the  same  class  of  disturbances,  i.e.,  an  initial  displace¬ 
ment  representing  the  response  to  a  static  load. 

Rew  and  Junkins'"  searched  for  optimal  su\te  and  control 
weighting  matrices  for  the  quadratic  performance  index  but  did  not 
consider  any  structural  variables.  Bodden  and  Junkins"  minimized 
some  unspecified  robustness  measure  while  placing  closed-loop 
eiaenvalues  in  a  desired  region.  When  the  structural  variables  were 
included,  a  sequential  approach  was  used.  McLaren  and  Slater*' 
employed  a  similar  homotopy  continuation  method  that  converged 
only  after  about  750  complete  analyses  for  one  case,  and  after  10  h 
of  CPU  time  on  a  Cray  Y-MP  supercomputer  for  another. 

Lust  and  Schmit'^  searched  for  structural  variables  and  linear 
output  feedback  gains  on  equal  footing  for  the  steady-state 
response  to  deterministic  harmonic  loading.  Special  attention  was 
paid  to  developing  explicit  approximations  to  the  highly  nonlinear 
implicit  response  functions.  Later,  Thomas  and  Schmit’’^  extended 
this  formulation  to  problems  with  noncollocated  actuators  and  sen¬ 
sors,  as  well  as  stability  constraints  on  the  complex  eigenvalues. 
Sepulveda  et  al.'^  extended  to  complex  eigenvalues  Canfield's 
Rayleigh  quotient  approximation  (RQA)  for  eigenvalues  of  con¬ 
servative  systems.'*  Based  on  the  complex  RQ.A.  Thomas  et  al. 
improved  the  approximations  for  damping  and  control  forces. 

Only  Belvin  and  Park'"  appear  to  have  taken  advantage  of  inde- 
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pendent  modal  space  control^  (IMSC)  to  treat  the  integrated  struc¬ 
tural/control  problem.  The  authors  did  not  claim  to  have  found  the 
integrated  optimum,  rather  they  used  IMSC  to  decouple  structural 
and  control  design  by  determining  the  primary  influence  of  the 
structural  design  on  the  control  objective. 

A  novel  approach  proposed  by  Messac  et  al.'^  combined  distur¬ 
bance  rejection  and  command  following  performance  errors  in  a 
quadratic  performance  norm.  One  of  the  final  designs  was  charac¬ 
terized  as  nonintuitive,  because  the  second  frequency  was  reduced 
by  half,  whereas  in  separable  design  structural  optimization  would 
normally  stiffen  the  second  mode,  increasing  its  frequency.'*  A 
comparison  of  designs  obtained  for  two  simple  truss  structures  was 
made  by  Rao^°  for  various  objective  functions. 


are  coefficient  matrices.  For  linear  feedback  of  the  controlled 
modes,  the  control  law  is 

/c(0=O^O«(0  =-Gr|c-'f/’ic  (6) 

where  the  modal  control  gain  matrices  G  and  H  are  constant  in  the 
steady-state  case.  To  implement  modal  feedback  control  the  modal 
displacements  and  velocities  must  be  extracted  from  the  system 
outputs.  Application  of  the  second  part  of  the  expansion  theorem^' 
provides  the  transformation  of  physical  to  modal  coordinates 

y\c{t)  =  (7a) 


Problem  Statement 

We  are  concerned  with  a  linear,  time-invariant,  distributed- 
parameter  system.  Following  spatial  discretization,  the  equations 
of  motion  are  given  by 

'  M{v;)q{t)  +  C{v^)q{t)  +  K{v,)q{t)=Du{v^,v^,t)  (1) 

where  q  is  the  configuration  vector,  M  the  mass  matrix,  C  the 
damping  matrix,  and  K  the  stiffness  matrix,  all  parameterized  by 
the  Wj  vector  Vj  of  structural  variables.  The  system  is  regulated  by  a 
vector  of  control  forces  «(Vj,  v^,,  i)  which  is  also  parameterized  by 
the  «;..-vector  v^.  of  control  design  variables,  and  transmitted  to  the 
structure  according  to  the  load  distribution  matrix  D.  The  struc¬ 
tural,  and  control  vectors  can  be  combined  into  the  complete 
design  vector  V  =  [v[  v].Y . 

Modal  Space  Control 

If  the  structure  is  subject  to  proportional  viscous  damping,  the 
uncontrolled  system  can  be  decoupled  by  means  of  the  classical 
modal  transformation.  By  the  expansion  theorem,*'  this  transfor¬ 
mation  is  given  by 


fic(0  =  (7b) 

in  which  the  q{t)  are  detenninistic  physical  coordinates.  State  esti¬ 
mation  is  not  considered  here;  however,  if  at  least  as  many  sensed 
outputs  are  available  as  controlled  modes,  modal  state  estimation 
is  not  an  obstacle.*  Meirovitch  et  al.**  successfully  demonstrated 
the  use  of  modal  filters  based  on  Eqs.  (7)  for  control  of  an  experi¬ 
mental  beam  by  means  of  IMSC. 

Although  they  are  not  controlled,  the  residual  modes  receive  the 
excitation 

/;;(r)  =  a>>«(/)  (8) 

The  closed-loop  system  is  described  by 

x=A,x  (9) 

where 


-G  -H 


g(.t)  =  (2) 

/  =  1 


where  <I>  is  the  classical  modal  matrix  and  T]  the  vector  of  modal 
coordinates.  The  eigenvectors  are  orthogonal  with  respect  to  the 
mass  matrix  and  can  be  normalized  so  as  to  satisfy  =  /, 

=  £2*,  where  £2*  is  the  diagonal  matrix  of  eigenvalues.  Syn¬ 
thesizing  a  control  law  is  more  convenient  when  the  left  side  of  Eq. 
(1)  is  transformed  to  diagonal  form.  To  this  end  we  insert  Eq.  (2) 
into  Eq.  (1),  premultiply  by  assume  that  damping  is  of  the  pro¬ 
portional  type,*'  and  obtain 

Ti  (r)  +  2^a  (V,)  fi (r)  +  £2'(v^)  p  (r)  =  v,,  r)=/  (t) 

(3) 

where  ^  is  the  diagonal  matrix  of  modal  damping  factors,  assumed 
to  be  constant  for  the  class  of  structures  under  consideration,  and 
J[t)  is  the  modal  force  vector. 

Typically,  not  all  n„  modes  need  be  controlled,  so  the  control 
design  problem  can  be  truncated.  To  this  end,  we  partition  the 
decoupled  equations  into  equations  for  the  controlled  modes 
and  ?}/(  equations  for  the  residual  modes.  Consistent  w'ith  this,  the 
modal  matrix,  the  modal  coordinate  vector,  and  the  modal  force 
vector  have  the  partitioned  forms  4>  =  [Oc  <!);;],  q  =  [q^ 
and  /(/)  =  [fcilY  fRitW,  respectively.  The  controlled  equa¬ 
tions  in  Eq.  (3)  can  be  transformed  to  state  space.  To  this  end,  we 
introduce  the  state  vector  x=  [qj^  q^ ]  .  Then,  the  state  equations 
for  the  controlled  modes  can  be  written  in  the  matrix  form 


where 


x  =Ax  +  Bu 


0  / 
-«c  -2Cci^c 


B  = 


0 


(4) 

(5) 


is  the  closed-loop  plant  matrix,  in  which 

G  =  £2c+G,  H  =  2I^c^c  +  H  (11) 

Independent  Modal  Space  Control 
The  general  idea  of  the  IMSC  method  is  that  the  control  force 
for  a  given  mode  depends  only  on  the  modal  displacement  and 
velocity  of  that  mode.  Hence,  the  independence  of  open-loop 
modal  equations  is  preserved  for  the  closed-loop  equations.  For 
self-adjoint  systems,  the  2/2^  ^  matrix  Riccati  equation  for 
linear  optimal  control  reduces  to  a  series  of  independent  2X2 
matrix  equations  for  each  of  the  controlled  modes.  It  follows  that, 
in  the  case  of  IMSC,  the  modal  gain  matrices  G  and  H  in  Eq.  (6) 
are  diagonal.  Equation  (6)  yields  the  modal  control  forces,  from 
which  we  obtain  the  actuator  forces 

u(t)={^cD)~'fciO  (12) 

When  the  number  of  controlled  modes  and  actuators  is  not  the 
same,  the  exact  inverse  (O^D)”'  in  Eq.  (12)  can  be  replaced  by  a 
pseudoinverse,  but  the  modal  forces  will  be  only  approximately 
independent,  depending  how  close  <I>jD  is  to  a  square  matrix. 

In  optimal  IMSC.  the  weighting  matrices  for  the  modal  space 
quadratic  performance  index  are  assumed  to  be  diagonal,  so  that 
the  performance  index  for  the  steady-state  case  can  be  expressed 
as  the  sum  of  independent  modal  performance  indices  in  the  form 

J  =  lj\x^Qx  +  fcRfc)  d/ 


(13) 
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As  a  result,  the  In^-  X  matrix  algebraic  Riccati  equation  re¬ 
duces  to  Hq  independent  2X2  matrix  algebraic  Riccati  equations. 
The  latter  have  the  analytic  solution 

+  '■r'  “ 

(14) 

/t,  =  7^'  +  -Sr  >  '■=  1.2,...,nc- 

One  choice  of  control  design  variables  is  the  modal  control 
force  weighting  factors  =  [i\r2  ■  ■  ■  because  they  determine 
the  modal  gains  for  IMSC  optimal  control  uniquely.  Alternatively, 
all  of  the  could  be  eliminated  by  using  Eqs.  (14)  to  solve  for 
each  /ir  in  terms  of  or  vice  versa.  In  the  former  case. 


ble  solution  exists.  Again,  parametric  variation  of  the  e,  generates  a 
set  of  Pareto  minima.  An  important  implication  for  the  problem 
under  consideration  is  that  any  single  objective  function  may  be 
minimized  whereas  the  others  are  treated  as  constraints.  The 
choice  of  criterion  is  at  our  discretion.  The  same  Pareto  optimal 
designs  should  be  generated  by  minimizing  either  structural 
w'eight  or  the  quadratic  performance  index,  for  example,  while 
constraining  the  other.  Of  course,  this  assumes  that  each  scalar 
optimization  finds  the  global  minimum  or  that  the  same  set  of  local 
minima  will  be  found  for  each  choice  of  objective  function.  The 
constraint  method  is  the  basis  for  seeking  Pareto  optimal  solutions 
of  this  research.  Thus,  the  focus  is  on  determining  efficient  meth¬ 
ods  to  solve  Eqs.  (18)  for  the  integrated  structural  and  control 
problem. 


+  2S,. 


(  2 


-  to. 


(15) 


whereas  in  the  latter  case. 


(16) 


Now,  either  the  g^  or  the  h,  can  be  taken  as  control  design  vari¬ 
ables. 


Pareto  Optimal  Solutions 

A  natural  approach  to  integrating  design  requirements  from  two 
or  more  disciplines  is  to  optimize  a  vector  of  cost  functions.  Ide¬ 
ally,  all  of  the  criteria  in  the  vector  would  be  optimized  simulta¬ 
neously,  Because  all  criteria  cannot  be  optimized  at  once,  an 
understanding  of  what  constitutes  a  vector  optimum  is  required. 
Following  Koski,-^  we  state  the  multiobjective  (multicriteria  or 
vector)  optimization  problem  as 

F*{v*)  =  min  F(v')  (17a) 

v£  n 

where  F:D.  — ^  is  a  s  ector  objective  function 

F(v)  =  Lf,(v)  F,(r)  ...  F„(r)J"  (17b) 

and  its  components  F;  :Q  —a  F,  /  =  1 . 2. . . . ,  m  are  the  criteria.  The 
design  variable  vector  v  belongs  to  the  feasible  set  Q  c  F  , 
defined  by  the  vectors  h  and  g  of  equality  and  inequality  con¬ 
straints,  respectively,  as 


n  =  {ve  F"|/i(v)  =  0,  g(v)  <0}  (17c) 

where  vector  inequalities  are  understood  to  apply  individually  to 
each  component  of  the  vector.  (Distinction  between  the  vectors  g 
and  h  for  constraint  functions  and  modal  gains  and  should  be 
clear  from  the  context  of  the  equations.)  In  general,  no  unique 
point  exists  that  optimizes  all  m  criteria  simultaneously.  Pareto 
defined  the  vector  optimum  v*  as  that  for  which  there  exists  no 
feasible  vector  r  that  would  decrease  some  criterion  without  simul¬ 
taneously  increasing  at  least  one  other  criterion. 

A  popular  technique  for  generating  Pareto  minima  is  the  con¬ 
straint  method.-^  In  this  method  the  vector  optimization  is  replaced 
by 

F't(v*)=  min  Fi^(v)  (18a) 

*  v€finn.(f) 

where 

n^Ce)  =  {v;F,(v)<£„  iaA-}  (18b) 

and 

ee  Fj=  {Le,e2...et_,ei,,  n(e)  =^0  (18c) 


Theoretical  Development 

Performance  Index  for  Integrated  Structural  and  Control  Design 

In  using  the  LQR  theory^  for  optimal  control,  one  determines  a 
control  law  that  minimizes  a  quadratic  performance  index  for  a 
given  positive  semidefinite  state  weighting  matrix  Q  and  a  positive 
definite  control  weighting  matrix  F.  The  choice  of  weighting 
matrices  is  otherwise  at  the  discretion  of  the  control  designer.  Each 
choice  produces  a  different  optimal  control.  In  practice,  the 
designer  alters  Q  and  F  to  balance  system  performance  and  control 
effort.  As  an  example,  one  particular  choice  might  result  in  satura¬ 
tion  of  an  actuator  in  some  simulation.  Then,  the  designer  might 
increase  the  weighting  factor(s)  associated  with  that  actuator.  In 
this  sense,  the  performance  index  is  not  used  to  compare  candidate 
control  laws.  Instead,  other  criteria  are  used  as  a  basis  for  compari¬ 
son  in  the  simulations.  The  dilemma  is  compounded  in  integrated 
structural  and  control  design,  in  which  performance  for  different 
plants  and  control  laws  must  be  compared.  Therefore,  the  question 
of  an  appropriate  performance  measure  for  controls  is  a  critical 
one.  We  propose  here  that  the  time  integral  of  the  total  system 
energy  be  used  as  the  single,  unique,  and  physically  meaningful 
performance  measure.  The  total  system  energy  consists  of  the 
kinetic  and  elastic  potential  energy  and  the  energy  expended  by  the 
control  system.  .As  mentioned  in  the  literature  review,  the  struc¬ 
tural  energy  has  been  used  frequently  to  define  a  unique  state 
weighting  matrix  Q.  Defining  the  control  weighting  matrix  F  is 
more  arbitrary.  Two  interpretations  of  “control  effort”  are  consid¬ 
ered  next. 

The  control  effort  in  the  performance  index  is  taken  to  represent 
the  energy  expended  by  the  control  system,  which  is  often 
assumed  to  be  proportional  to  the  square  of  the  control  inputs. 
When  the  input  represents  an  actuator  force,  the  energy  to  actuate 
certain  mechanical  systems  might  be  proportional  to  the  square  of 
each  actuator’s  force.  In  this  case,  we  expect  the  energy  expended 
by  each  actuator  to  be  independent  of  the  others,  in  which  case  F  is 
diagonal.  Assuming  that  the  proportionality  factors  of  each  actua¬ 
tor  are  known  a  priori,  we  consider  a  performance  measure  for^the 
system  defined  by  Eq.  ( 1 )  in  the  form 


where  F  is  a  unique  weighting  matrix  of  factors  for  the  control 
energy. 

Before  proceeding  any  further,  let  us  consider  alternative  mea¬ 
sures  for  the  control  effort.  The  weighting  matrix  F  =  was 

proposed  by  Venkayya  and  Tischler*  and  employed  by  Rao“°  and 
Belvin  and  Park.'*  This  weighting  matrix  provides  a  quasistatic 
measure  of  the  energy  imparted  to  the  structure  by  the  control  sys¬ 
tem.  It  can  be  derived  by  assuming  that  the  forces  are  applied  qua- 
sistatically  to  the  structure  by  the  control  Du{t).  The  elastic  energy 
imparted  to  the  structure  has  the  expression 

{\/2)[Du{t)\^q{t)  =  (1/2)  [Du{t)VK'Du{t) 


One  criterion  is  taken  as  a  scalar  objective  function  whereas  the  ^ 

remaining  are  con.strained  by  a  set  of  constants  £,  such  that  a  feasi-  =  (^/2.)u{t)Ru(t)  (20) 
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where  q{t)  represents  the  static  displacement  due  to  the  quasi¬ 
static  force  Du{f).  However,  one  could  as  well  consider  the  energy 
imparted  to  the  structure  by  the  control  force  vector  Du{t)  applied 
impulsively.  Both  are  approximate  measures  of  the  actual  energy 
used  by  the  control  system  to  do  work  on  the  structure. 

Next,  we  consider  the  actual  work  done  by  the  actuator  forces, 
or 

ly  =  f’/F(/)dr  =  f'fiV(0df  =  (21) 

Jo  Jo 


conservative  plant  the  useful  work  is  simply  the  change  in  internal 
energy  of  the  plant.  Therefore,  when  comparing  two  systems,  it  is 
sufficient  to  use  the  unique  state  weighting  matrix  that  produces 
the  total  structural  energy  and  let  =  0  when  comparing  two  sys¬ 
tems.  Of  course,  the  second  choice  ignores  control  spillover,  which 
must  be  small  for  any  feasible  system.  We  also  assumed  the  con¬ 
trol  gains  were  such  that  the  actuators  did  not  reach  saturation. 
Next,  we  derive  the  time-invariant  form  of  the  performance  index 
and  formulate  constraints  for  the  actuator  forces. 


where 


=  [  Tic/c(^)d'.  =  [  Ti«/«(i:)dr  (22) 

Jo  Jo 

are  the  work  done  by  the.  controlled  modes  and  by  the  residual 
modes,  respectively.  At  this  point,  we  define  the  useful  work  as  the 
work  done  by  the  actuator  forces  on  the  controlled  modes,  i.e., 
and  note  that  for  IMSC  we  can  evaluate  analytically.  Because 
thl5  modal  forces  are  linear  in  the  state  variables,  we  assume  that 
the  integral  has  a  quadratic  form,  or 

f'Tic(^)/c(^)dt  =  a:^(t)Fx(/)  (23) 

Jo 

in  which  x  =  fie is  the  modal  state  vector.  Differentiating 
both  sides  of  Eq.  (23)  with  respect  to  time,  we  have 

tic  (0  /c  (0  =  x\t)  Rx{t)+  x\t)  Rx{t)  (24) 

Substituting  X  =  tic  2nd  (6)  into  the  left  side  of  Eq.  (24) 
and  Eq.  (9)  into  the  right  side,  we  obtain 


0  -(1/2)G 

-(1/2)G  -H 


X  =  x^  [a]r  +  RA  ]x 


(25) 


Equation  (25)  will  be  satisfied  for  an  arbitrary  state  only  if  the 
matrices  in  brackets  are  equal.  Because  for  IMSC  the  various  sub¬ 
matrices  are  diagonal,  insertion  of  Eq.  (10)  into  Eq.  (25)  yields 


Independent  Model  Space  Control  Time-Invariant  System  Energy 
We  wish  to  minimize  the  time  integral  of  the  system  energy,  Eq. 
(19),  as  a  function  of  the  vector  of  structural  parameters  Vj  deter¬ 
mining  the  natural  frequencies  and  the  vector  of  control  parame¬ 
ters  determining  the  control  gains.  Because  the  plant  and  control 
are  taken  to  be  time  invariant,  the  performance  index  can  be 
expressed  without  regard  to  time,  assuming  that  the  closed-loop 
system  is  stable.  The  perfonnance  index.  Eq.  (19).  can  be  written 
in  the  partitioned  modal  state  space  form 


jric(0") 

Uc(0") 


"  Pu  Pm  jfic(0")| 

Pm  P22  Uc(0^)J 


(28) 


The  positive  definite  symmetric  coefficient  matrix  P  is  a  solution 
to  the  Lyapunov  equation 

Q,+  aIp+PA,  =  0  (29) 


for 


Qc 


Q.-  0 
0  I 


+  IG,  H^]‘R[G^ 


(30) 


w'here  G^,  and  relate  the  actuator  forces  to  modal  coordinates 
according  to 


M(0  =  -G^qc-.W,iTic 


(31) 


gh^'h-g 

0 


h-'h 


(26) 


where  G  and  H  are  defined  by  Eqs.  (11).  When  no  natural  damp¬ 
ing  is  present  H  =  H  and  Eq.  (26)  reduces  to 


0 

/ 


(27) 


The  use  of  IMSC  in  formulating  Eqs.  (28-30)  has  significant 
implications.  First,  the  actuator  gains  can  be  synthesized  from  the 
modal  control  gains  by  means  of  Eq.  (6)  or 

G„  =  ((I>^D)"'g,  //^  =  (O^D)"’//  (32) 

As  a  result,  even  if  R  is  diagonal  in  Eq.  (30),  the  modes  recouple 
Eq.  (29)  via  Eqs.  (32).  However,  choosing  /?  =  0  becomes  quite 
attractive  when  using  IMSC.  Then  becomes  diagonal,  thus 
decoupling  Eq.  (29)  and  permitting  an  analytical  determination  of 
the  diagonal  partition  matrices  of  P.  Indeed,  inserting 


which  is  precisely  equal  to  the  modal  state  weighting  matrix  Q  ren¬ 
dering  {\l2)x^Qx  the  sum  of  kinetic  and  strain  energy.  This  result 
is  to  be  expected,  because  for  a  conservative  structural  system  the 
change  in  energy  from  the  initial  time  to  a  quiescent  state  at  the 
final  time  must  be  dissipated  entirely  by  the  control  system. 
Hence,  if  we  define  control  effort  as  the  useful  work  done  by  the 
control  system,  it  is  sufficient  to  merely  include  the  state  weighting 
matrix  in  the  performance  index.  It  is  important  to  stress  here  that 
Eq.  (19)  with  P  =  0  represents  a  performance  index  for  comparing 
two  systems  and  not  a  performance  index  used  to  derive  LQR  con¬ 
trol  gains. 

In  summary,  in  seeking  an  appropriate  and  physically  meaning¬ 
ful  control  performance  measure  we  arrived  at  two  alternatives. 
One  employs  the  unique  matrix  R  in  Eq.  (19)  yielding  the  energy 
expended  by  the  actuators  over  time,  assuming  that  such  a  charac¬ 
terization  of  the  control  system  is  available.  In  the  absence  of 
knowledge  about  R,  the  second  alternative  is  to  consider  the  useful 
work  done  by  the  control  system  over  a  given  time  interx'al.  For  a 


11  '12 
P 

12  ^^22 


(33) 


and  Eqs.  (10),  (11),  and  (30)  in  to  Eq.  (29),  we  obtain 


co‘ (2^,03,.  + /;,-)  2a)" -tg. 


2(co,'-(-g,.)  2(2C,co,-F/t,) 


(34a) 


0), 


P12,  = 


2(co"-tg,) 


(34b) 


2co-  +  g, 


2(a)--Hg,)(2C,a),  +  /7,) 


/  =  1, 2, . . . ,  (34c) 
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A  critical  assumption  must  be  examined  for  properly  minimiz¬ 
ing  the  integrated  objective  function.  Eq.  (28).  An  arbitrary  exter¬ 
nal  disturbance  for  a  fixed  plant  is  often  assumed  to  induce  an  arbi¬ 
trary  initial  state  of  unknown  direction  but  constant  magnitude. 
Minimizing  Eq.  (28)  for  this  case  is  then  shown  to  be  equivalent  to 

min  J  <=>  min  [tr(P)]  (35) 


was  derived  by  differentiating  Eq.  (34c).  The  structural  frequency 
derivative  in  Eq.  (40)  is  well  known.'  Note  that  its  coefficient  on 
the  right  side  of  Eq.  (40)  is  positive. 

— ?>0  (41) 

dco" 


In  earlier  approaches  to  the  integrated  problem,  this  assumption  has 
often  been  made  in  a  casual  manner.  However,  if  design  of  the  plant 
itself  is  under  consideration,  the  initial  state  depends  on  the  struc¬ 
tural  properties.  Indeed,  an  arbitrary  external  disturbance  of  given 
magnitude  produces  a  corresponding  initial  state  depending  on  the 
as  yet  undetermined  stiffness  and  mass  of  the  structure.  Therefore, 
the  sensitivity  of  the  initial  state  d.r:g/di\.  cannot  be  ignored.  Cal¬ 
culating  dxg/dv^.  involves  the  assumption  of  a  certain  disturbance, 
e.g.,  an  impulsive  force  of  given  magnitude  or  at  least  assumed  sta¬ 
tistics  governing  the  nature  of  the  disturbance.  The  same  assump¬ 
tions  apply  to  calculating  actuator  forces,  as  well.  For  a  given  dis¬ 
turbance,  the  derivative  of  the  performance  index,  Eq.  (28),  with 
respect  to  structural  and  control  variables  is 


dv. 


T  dP 


T 


(36) 


Sensitivity  of  the  first  term  requires  differentiation  of  Eq.  (29).  The 
second  term  in  Eq.  (36)  involves  structural  variables  only.  Both 
terms  are  included  in  the  following  derivation  for  the  IMSC  perfor¬ 
mance  index. 

Consider  an  impulsive  disturbance  represented  by  the  force  vec¬ 
tor  where  Q  is  the  intensity  vector  and  5(f)  is  the  Dirac 

delta  function.  Solution  of  the  modal  counterpart  to  Eq.  (1)  with 
the  right  side  replaced  by  the  impulsive  force  25(f)  results  in  the 
initial  modal  velocities 


fi(0")  =  <d’’(v,)2  (37) 

By  substitution  of  Eq.  (37)  into  Eq.  (28),  the  performance  index 
can  be  simplified  to 

J  =  fi’'(0")/>3jfi''(0")  = 


so  that,  ignoring  the  modal  sensitivity  appearing  in  Eq.  (39),  we 
conclude  that  a  reduction  of  the  performance  index  implies  a  reduc¬ 
tion  of  structural  frequencies,  i.e.,  a  softening  of  the  structure.  This 
result  is  in  opposition  to  the  heuristic  conclusion  reached  by  Belvin 
and  Park,'*  but  agrees  with  the  observation  of  Messac  et  al.'^ 
Differentiation  of  Eq.  (38)  with  respect  to  the  control  design 
variables  leads  to 


dp2 

d\’ 


(42) 


The  partial  derivative  with  respect  to  a  control  design  variable  is 
the  partial  derivative  with  respect  to  the  IMSC  modal  rate  gain, 

=  hi..  Hence,  by  differentiating  Eq.  (34c)  and  using  Eq.  (16)  for 
optimal  IMSC  gains,  the  partial  derivative  on  the  right  side  of  Eq. 
(42)  can  be  shown  to  have  the  form 


^P22.  -CO~  C0^  +  g;/2 

2  (cof  +  gf  JTThfTaf  +  g,.) 


(43) 


and  note  that  9^,, /0/!,.  <  0  so  that  decreasing  the  performance 
index  implies  increasing  the  modal  gains,  as  expected.  However, 
the  modal  gains  will  be  prevented  from  becoming  infinitely  large 
by  constraints  on  the  actuator  forces  posed  in  the  next  section. 


.Actuator  Force  Constraints 

Actuator  forces  are  constrained  to  be  smaller  than  some  maxi¬ 
mum  magnitude  for  all  time.  To  create  a  finite  number  of  con¬ 
straints  for  the  integrated  optimization,  the  actuator  inputs  (forces) 
can  be  evaluated  at  only  a  few  peak  times.  Peak  actuator  forces 
occur  either  at  r  =  0  or  when  the  derivative  of  the  force  with  respect 
to  time  vanishes.  In  the  first  case,  the  peak  inputs  for  each  actuator 
are  easily  determined  from  «(0'^)  =  -(jXq  .  Inserting  Eqs.  (6)  and 
(37)  into  Eq.  (12),  the  initial  IMSC  actuator  forces  corresponding 
to  an  initial  velocity  caused  by  an  impulse  of  intensity  Q  are 


^  ^  ■  I  A  7 

=  <& 

/  =  I 


(38) 


^  J  •. 

where  ((),  is  the  ith  component  of  the  vector  <t>  =  0(-2  -  so  that  the 
control  objective  represented  by  Eq.  (38)  a  simple  weighted  sum  of 
the  diagonal  elements  of  ‘’c )•  It  should  be  pointed  out  here 

that  an  initial  condition  due  to  an  impulse  is  more  realistic  than  a 
static  preload  and  applies  to  unconstrained  as  well  as  to  con¬ 
strained  structures.  Sensitivity  of  the  integrated  IMSC  perfor¬ 
mance  index  for  an  impulsive  disturbance,  Eq.  (38),  is  considered 
in  the  next  section. 


Independent  Modal  Space  Control  Performance  Index  Derivatives 
Differentiation  of  Eq.  (38)  with  respect  to  the  structural  vari¬ 
ables  leads  to 


K 


9  A 

dv„ 


(39) 


u(0+)  =  -(a)cD)-‘//Oce  (44) 

Subsequent  peak  inputs  first  require  determination  of  the  peak 
times  from  the  transcendental  equation  resulting  from  differentiat¬ 
ing  the  scalar  form  of  Eq.  (12) 

=  -G,A^.x{l)  ,  V  G  .N,,  (45) 

where  Gy  represents  the  jth  row  of  the  gain  matrix  G.  Solution  of 
Eq.  (45)  generates  a  set  of  Pj  peak  times  Ty  ={?[,...  ,  t^y )  for  each 
of  the  /(„  inputs.  For  IMSC  the  modal  states  are  known  explicitly  in 
terms  of  the  design  variables  and  time."'' 

Actuator  Force  Derivatives 

Consider  first  the  actuators  forces  at  time  t  =  due  to  an 
impulse.  The  initial  state  vector 


The  partial  derivative  term  3^,  /0Vjy  in  Eq.  (39)  requires  natural 
eigenvector  derivatives.  An  explicit  expression  for  the  partial 
derivative  term 


^Pii.  hi  dco~ 

4  ( cof  -F  g,  )~J4  +  hf/(o^ 


(40) 


does  not  vary  with  time  or  control  variables,  only  with  structural 
variables.  Therefore,  we  divide  the  possible  derivatives  into  those 
that  depend  on  structural  variables  and  those  that  depend  on 
control  variables  Differentiating  Eq.  (44)  with  respect  to  struc- 
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tural  variables  first  and  noting  that  H  does  not  depend  explicitly  on 
Vjj ,  we  have 


9u(0^) 


d(^cQ)  d0r  T  ->  /fj  n 


(47) 


k=l,...,n. 


where  the  modal  gradients  were  first  encountered  in 

Eq.  (39).  Differentiation  of  Eq.  (44)  with  respect  to  control  vari¬ 
ables  yields 


ov. 


A' =  1, ....  n. 


(48) 


Table  1  Independent  modal  space  control  gains 
for  simply  supported  beam 


Case 

Mode  I 

P, 

Si 

A, 

C/ 

Initial 

1 

1.0 

0.4987 

1.4133 

0.0714 

IMSC 

2 

I.O 

0.4999 

1.4142 

0.0179 

3 

I.O 

0.5000 

1.4142 

0.0080 

Optimized 

I 

0.1769 

2.7874 

3.3510 

0.1674 

IMSC 

2 

0.1391 

3.5909 

3.7910 

0.0480 

3 

0.1639 

3.0500 

3.4930 

0.0197 

Optimized 

1 

0.1955 

2.517 

3.186 

0.1782 

Structure 

2 

0.0989 

5.045 

4.494 

0.0614 

and  IMSC 

3 

0.1505 

3.322 

3.646 

0.0199 

The  modal  rate  gains  along  the  diagonal  H  are  themselves  the  con¬ 
trol  design  variables,  so  that  dH is  a  matrix  of  all  zeros,  ex- 
c^t  for  entry  in  the  Ath  row  and  column  which  is  equal  to  one.  Ac¬ 
tuator  force  derivatives  at  subsequent  peak  times  can  also  be 
derived, but  were  not  considered  in  the  examples  to  follow. 

Example:  Simply  Supported  Beam  with 
Three  Actuators 

The  performance  index,  Eq.  (38),  was  minimized  for  the  pur¬ 
pose  of  designing  the  IMSC  modal  gains  for  a  uniform  beam  with 
three  actuators,  described  in  Ref.  3  (Fig.  1).  The  material  proper¬ 
ties  are  given  in  unspecified  consistent  units  such  that  EII(mL^)  = 
1.  The  disturbance  was  a  unit  impulse  at  0.43Z,.  The  force  actuators 
were  located  0.15L,  0.55L,  and  0.73L  of  the  span.  Coupled  control 
to  damp  out  the  first  eight  modes  produced  a  maximum  actuator 
force  of  magnitude  3.8  units. ^  The  first  12  modes  were  used  to 
simulate  the  closed-loop  response  of  the  beam  to  the  unit  impulse. 
Time  histories  of  the  midspan  deflection  of  the  beam  and  the  cor- 


Fig.  3  Pareto  optima  for  simply  supported  beam. 


Impulse 


1  T 


I  F(t)  =  actuator  forces 
Fig.  1  Simply  supported  beam  with  three  actuators. 


Fig.  2  Time  histories  of  beam  with  three  actuators. 


responding  actuator  forces  appear  in  Fig.  2  for  three  cases  using 
IMSC,  described  next. 

IMSC  was  used  to  control  only  the  first  three  modes.  For  unit 
values  of  the  modal  control  weighting  factors  in  Eq.  (13),  actua¬ 
tor  forces  were  lower  and  their  time  history  smoother  than  for  cou¬ 
pled  control,  but  at  the  expense  of  decreased  damping.  Next,  IMSC 
modal  gains  were  treated  as  design  variables  and  Eq.  (28)  was 
minimized  subject  to  constraints  on  each  actuator  force.  The  upper 
limit  on  the  actuator  forces,  3.8  units,  was  taken  as  the  maximum 
value  from  the  original  example.  Only  a  single  call  to  an  IMSL 
nonlinear  optimization  subroutine  was  needed.  Results  in  Fig.  2b 
demonstrate  higher  damping  and  a  14%  reduction  in  the  maximum 
deflection  without  saturating  the  actuators.  The  modal  control 
weighting  factors  and  IMSC  modal  gains  for  each  case  are  given  in 
Table  1.  This  case  represents  the  best  performance  possible  by 
optimizing  control  gains  alone  without  design  structural  variables. 
In  other  words,  it  is  the  optimal  control  for  a  fixed  structure. 

Next,  structural  design  variables  were  introduced.  To  establish 
the  relationship  between  bending  and  inertia  properties  (i.e.,  El  and 
niL),  a  particular  cross-sectional  shape  was  specified.  A  circular 
tube  of  outer  radius  R  =  Ll\00  and  thickness  r  =  R/IO  was  consid¬ 
ered.  TTiickness  was  selected  as  the  structural  design  variable  while 
the  outside  radius  remained  fixed.  The  minimum  thickness  allowed 
was  r^in  =  R/4Q,  whereas  the  maximum  thickness  was  f^ax  = 

To  approximate  the  continuous  distribution  of  material  along  the 
span,  an  Euler-Bemoulli  beam  was  modeled  using  100  finite  ele¬ 
ments.  Each  structural  design  variable  controlled  the  thickness  of  a 
group  of  elements  located  symmetrically  about  the  midspan.  The 
performance  index,  Eq.  (38),  w'as  minimized  w'hile  requiring  no  in¬ 
crease  in  total  mass  and  limiting  the  actuator  forces.  Also,  the  fun¬ 
damental  natural  frequency  was  constrained  to  be  at  least  80%  of 
its  initial  value.  The  integrated  optimization  problem  may  be  stated 
as  minimizing  the  system  energy 

min  7(v,.,  v,) 

V  .  V 


(49a) 
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Fig.  4  Optimum  beam  profiles  for  10  structural  design  variables. 


Fig.  5  Optimum  beam  profile  for  50  structural  design  variables. 


Structural  design  variables  v,  (tube  thickness)  and  control  design 
variables  v^.  (modal  gains)  were  determined  simultaneously.  The 
time  history  in  Fig.  2c  indicates  a  29%  reduction  in  ma.ximum 
amplitude  without  loss  of  damping  or  exceeding  the  actuator  limit. 
This  third  case  demonstrates  the  benefit  of  multiobjective  optimi¬ 
zation.  Performance  improvement  over  Fig.  2b  was  achieved  by 
simultaneously  determining  control  design  variables  and  structural 
design  variables. 

Pareto  optimal  solutions  were  generated  by  parametrically  vary¬ 
ing  the  mass  constraint.  The  resulting  Pareto  optimal  sets  in  Fig.  3 
show  that  not  many  structural  design  variables  (SDV)  were  needed 
to  approximate  the  Pareto  optimal  objective  function  values.  Re¬ 
sults  for  more  than  10  structural  design  variables  were  not  shown 
because  they  were  indistinguishable  from  the  curve  for  10  vari¬ 
ables.  The  same  Pareto  curves  were  generated  by  minimizing  the 
mass  subject  to  a  constrained  performance  index.  The  resulting 
beam  profiles  in  Fig.  4  for  minimizing  the  performance  index  with 
10  structural  design  variables  were  identical  to  those  found  by 
minimizing  mass.  Of  course,  the  true  optimum  structural  design 
entails  a  continuously  varying  cross  section,  simulated  in  Fig.  5  by 
using  50  structural  design  variables  for  mass  constrained  to  unity 
(a  refinement  of  Fig.  4b). 

The  effect  of  the  eigenvector  derivative  terms  appearing  in  Eq. 
(42)  is  demonstrated  by  the  iteration  histories  given  in  Fig.  6.  The 
curve  labeled  “SQP”  represents  a  solution  of  the  integrated  optimi¬ 
zation  problem  using  sequential  quadratic  programming  to  solve 
Eqs.  (49)  directly.  The  curve  labeled  “Approx.  Concepts"  is  a 
more  efficient  solution  in  which  all  structural  and  control  functions 
are  approximated  by  first-  or  second-order  Taylor  series.^'*  The 
resulting  approximate  subproblem  was  solved  iteratively  until  it 
converged  to  the  exact  solution  found  by  SQP.  The  third  curve, 
labeled  “No  Eigenvector,"  was  an  attempt  to  solve  the  integrated 
problem  while  ignoring  the  first  term  in  the  summation  of  Eq.  (42). 
Without  the  eigenvector  sensitivity,  the  solution  converged  to  a 
distinctly  different  nonoptimal  design  with  higher  system  energy. 
This  alternate  design  was  achieved  only  when  using  approxima¬ 
tion  concepts.  The  direct  SQP  algorithm  failed  to  converge  at  all 
when  eigenvector  derivatives  were  ignored. 


Fig.  6  Iteration  histories  for  mass  =  1. 


subject  to  behavior  constraints  on  nondimensional  mass,  funda¬ 


mental  frequency,  and  actuator  constraints. 

m(v,)<l  (49b) 

a)'(v,)  >0.8(co])o  (49c) 

My(v„  v,,0  <3.8,  y=  1,2,3  (49d) 

as  well  as  side  constraints 

R/40<v,i,<R/5,  k=l,2,...,l0  (49e) 

0<va,  lc=\,2,3  (490 


Conclusion 

The  integrated  structural  and  control  design  problem  was  posed 
as  a  multiobjective  optimization  problem.  The  total  energy  inte¬ 
grated  over  time  was  identified  as  the  performance  index  for  con¬ 
stant  feedback  gain  closed-loop  control.  Independent  modal  space 
control  enabled  a  closed-form  solution  of  the  system  energy  as  a 
function  of  open-loop  frequencies,  modal  damping,  and  modal 
gains.  Peak  actuator  forces  were  constrained  to  be  within  pre¬ 
scribed  limits.  Pareto  optima  for  the  multiobjective  optimization 
problem  were  generated  using  the  constraint  method.  The  resulting 
scalar  optimization  problem  was  efficiently  solved  by  forming  an 
explicit  approximate  subproblem.  Results  revealed  that  eigenvec¬ 
tor  sensitivity  was  important  to  converging  to  the  integrated  opti¬ 
mum.  More  importantly,  an  integrated  design  approach  clearly 
improved  the  performance  of  the  closed-loop  system  without  an 
increase  in  structural  mass.  Once  the  superiority  of  the  integrated 
approach  was  established,  Pareto  optimal  curves  were  developed 
to  illustrate  the  extent  to  which  structtirai  mass  changes  could 
affect  the  closed-loop  performance  index. 

Several  important  insights  were  gained  from  this  research.  First, 
a  unique  and  physically  meaningful  performance  index  such  as  the 
total  system  energy  lends  it.self  to  a  fair  comparison  among  candi¬ 
date  structural  and  control  systems.  Furthermore,  IMSC  offered 
the  advantage  not  only  of  computational  efficiency,  but  more 
importantly,  it  permitted  the  derivation  of  explicit  expressions  for 
the  performance  index  and  its  derivatives.  We  infer  from  these 
derivatives  that  the  total  energy  in  the  system  imparted  by  an 
impulsive  disturbance  is  reduced  by  softening  the  structure,  i.e.,  by 
reducing  its  natural  frequencies.  Moreover,  the  eigenvector  deriva¬ 
tives  account  for  the  effect  structural  changes  have  on  the  initial 
condition  resulting  from  fixed  external  disturbances.  The  open- 
loop  modes  also  determine  the  effectiveness  of  actuators.  Thus, 
structural  changes  may  help  satisfy  constraints  on  actuator  forces. 
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Explicit  actuator  force  constraints  were  incorporated  as  the  physi¬ 
cal  means  by  which  control  gains  were  bounded.  This  approach  of 
minimizing  the  system  energy  and  structural  mass  subject  to  con¬ 
straints  on  frequencies  and  actuator  forces  made  integrated  struc¬ 
tural  and  control  design  a  tractable  problem.  The  use  of  IMSC  and 
the  selection  of  modal  gains  as  the  control  design  variables  were 
the  key  to  producing  computational  solutions  for  Pareto  optimal 
designs,  thus  demonstrating  the  tradeoff  between  structural  and 
control  objectives. 
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Abstract 

The  equations  of  motion  for  a  nonuniform,  aniso¬ 
tropic  thin-walled  beam  are  derived  and  applied  to 
the  study  of  vibration  and  static  aeroelastic  instabil¬ 
ity  of  slender  tapered  aircraft  wings  made  of  advanced 
composite  materials.  Numerical  results  illustrate  the 
effects  of  anisotropy.^transverse  shear  flexibility,  pri- 
meiry  and  secondary  warping,  as  well  as  of  wing  taper 
ratio,  and  the  implications  of  these  effects  on  the  vi¬ 
brational  and  divergence  instability  characteristics  are 
discussed. 

Introduction 

The  demands  on  advanced  flight  vehicles  to 
achieve  the  highest  performance  possible  over  the 
broadest  range  of  missions  has  stimulated  a  great  deal 
of  research  lately.  Improvement  of  flight  vehicle  struc¬ 
tural  efficiency  can  be  achieved  by  reducing  weight, 
controling  wing  deformations,  as  well  as  the  aeroelas¬ 
tic  load  distribution,  increasing  the  flutter  speed,  in¬ 
creasing  control  effectiveness,  etc.  Through  the  use  of 
advanced  composite  materials,  the  designer  acquires  a 
great  deal  of  freedom  in  carrying  out  structural  and 
aeroelastic  tailoring.  This  amounts  to  taking  advan¬ 
tage  of  the  light  weight  and  the  cross-couplings  of  com¬ 
posite  materials  to  optimize  the  response  characteris¬ 
tics. 

To  permit  results  useful  in  the  design  of  advanced 
flight  vehicles,  the  structured  models  must  be  as  com¬ 
prehensive  as  possible,  which  implies  that  they  must 
include  all  the  essential  characteristics  of  real  eurcraft 
structures.  To  meet  this  requirement,  in  a  number  of 
recent  studies^"®,  a  refined  wing  structural  model  of 
advanced  flight  vehicles,  including  many  effects  unex¬ 
plored  heretofore,  was  developed  and  the  implications 
of  the  newly  introduced  effects  were  revealed. 

One  factor  that  has  not  received  sufficient  atten¬ 
tion  to  date  is  the  nonuniformity  of  the  wing  cross 
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section.  It  seems  that  for  composite  aircraft  wings 
modeled  as  solid  beams,  the  only  results  for  nonuni¬ 
form  wings  were  presented  in  Refs.  3  and  6.  However, 
there  appears  to  be  no  investigation  of  wing  structures 
modeled  as  nonuniform  thin-walled  beams.  It  is  one  of 
the  goals  of  this  paper  to  develop  a  tapered  thin-wdled 
beam  model  for  slender  swept  wings. 

Due  to  renewed  interest  in  flight  vehicles  operat¬ 
ing  at  supersonic  and  hypersonic  speeds,  a  biconvex 
cross-section  profile  of  the  wing  will  be  considered. 
In  addition  to  cross-sectional  nonuniformity,  the  thin- 
waUed  wing  model  will  also  incorporate  anisotropy  and 
heterogeneity  of  the  composite  structure,  transverse 
shear  flexibility  of  the  constituent  materials  and  pri¬ 
mary  and  secondary  warping  effects.  Based  on  the 
anisotropy  of  the  material,  the  ply  angles  inducing  the 
most  favorable  elastic  cross-couplings  from  the  point 
of  view  of  aeroelastic  2uid  dynamic  response  character¬ 
istics  wiU  be  implemented. 

The  equations  of  motion  and  associated  bound¬ 
ary  conditions  for  wing  structure  models  as  described 
above,  amd  permitting  study  the  static  aeroelastic  be¬ 
havior  2Uid  vibrational  response,  Me  derived  by  means 
of  the  extended  Hamilton’s  principle.  Upon  solving 
the  related  eigenvalue  problems,  the  influence  of  the 
cross-sectional  nonuniformity  and  of  other  nonclassi- 
cal  features  is  investigated. 

Basic  Assumptions 

Under  consideration  is  a  structural  model  consist¬ 
ing  of  a  thin-walled  beam  and  intended  to  simulate  the 
lifting  surface  of  advanced  flight  vehicles.  The  cross- 
sectionad  contour  has  a  biconvex  profile,  with  geomet¬ 
rical  characteristics  varying  linearly  along  the  wing 
span,  (see  Fig.  1).  The  model  incorporates  the  follow¬ 
ing  nonclassical  features  (some  of  them  used  in  Ref.  5): 

i)  Anisotropy  of  constituent  materials, 

ii)  Transverse  shear  flexibility, 

iii)  Primary  and  secondary  warping  effects  and 

iv)  A  tapered  wing  surface  with  geometrically  similar 

cross  sections  at  aU  spanwise  stations. 

The  latter  requires  the  following  linear  distribution 
along  the  span  of  the  chord  c{t)),  height  6(rj)  of  the  mid¬ 
line  cross-section  profile  and  of  the  offset  e{r])  between 
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the  aerodynamic  and  reference  axes 

(civ)]  fcii] 

{  biv)\  =  [l-Vil-<T)]{  bR  \  (1) 

[civ)}  [ch) 

where  <t~=  ct/cr  defines  the  taper  ratio,  j;  =  r/£ 
is  the  dimensionless  spanwise  coordinate,  where  L 
denotes  the  wing  semi-span,  and  the  subscripts  R  and 
T  refer  to  the  root  and  tip  wing  sections.  Moreover, 

iZ(7?)  =  [!-»,  (l-<T)]iZfl  (2) 

in  which  i?(j7)  is  the  radius  of  curvature  of  the  circular 
arc  associated  with  the  midline  contour  at  section  v, 
Rr  is  the  radius  of  curvature  of  the  root  section  of 
the  wing.  Note  that  Qpints  on  the  beam  cross  sections 
are  identified  by  the  global  coordinates  x,y,  while  z 
is  the  spanwise  coordinate.  The  assumption  of  in- 
plane  nondeformability  of  the  beam  cross  section  is 
also  adopted. 

Due  to  the  geometrical  similarity  of  the  cross 
sections,  the  angles  0  between  the  axis  y  and  the  radii 
R  taken  at  points  cr/2,  c/2  and  ct/2  are  equal,  i.e. 
Qr  =  Qt  =  0  =  ©o- 


functions  not  only  of  the  circumferential  coordinate 
s  but  also  of  the  spanwise  coordinate  z.  As  a  result, 
free  warping  is  precluded  for  the  problem  at  hand. 

In  accordance  with  i)-iii),  and  in  order  to  reduce 
the  three-dimensional  elasticity  theory  of  beams  to  an 
equivalent  one-dimensional  one,  the  components  of  the 
displacement  vector  are  expressed  as 


u  (z,  y,z,t)  =Ug  (r,  t)  -  y  (s,  z)  4>  (z,  t) 
V  (z,  y,  z,  t)  =Vo  iz,  t)  +  xis,z)<i>iz,t) 


wix,y,z,t)  =Woiz,t)  +er  (z,t) 


y(s,z)  -  n 


a:(s,z)  +  n^ 


+  Sy  (^lO 

-  (».  2^)  +  na  (s.  z)] 


(5a) 

(56) 

dx 

ds 


(5c) 


in  which  («,  z)  and  no  (s,  z)  are  related  to  the  pri¬ 
mary  and  secondary  warping  functions,  respectively, 
and  n  denoting  the  coordinate  in  the  thickness  direc¬ 
tion.  In  addition. 


a  (s,  z)  =  -y  (s,  z)  dy/ds  -  z  (s,  z)  dy/ds  (6a) 


and 


Kinematical  Equations 

In  view  of  features  iii)  and  iv),  the  primary  warp¬ 
ing  function  becomes  a  function  of  both  s  and  q, 
i.e., 

F^  =  F^J  (s,  v)=  f  [^n  (s)  -  ds  (3a) 
Jo 

where 

denotes  the  torsiond  function,  in  which  Adv)  = 
/  r„  (s,  T])  ds  is  the  cross-sectional  area  of  the  beam 
bounded  by  the  midline  contour  at  rj,  sir])  denotes  the 
arc-length  measured  along  the  circumferential  coordi¬ 
nate  (whose  origin  is  arbitrarily  but  conveniently  cho¬ 
sen),  s  is  a  dummy  coordinate  associated  with  the  s- 
coordinate,  §i-)ds  denotes  the  integral  along  the  closed 

midline  contour,  r„  (s,  »?)  =  *(«,»?)  ^  -  y  («,  q)  ^  de- 

ds  ds 

notes  a  geometric  quantity  (see  Ref.  5)  and  /^(q)  = 
/  ds  iv)  is  the  circumferential  contour  length  at  q.  Af¬ 
ter  lengthy  calculations  one  obtains 

i>iv)  =  R  iv) 

Ac  iv)  =  R^  iv)  (200  -  sin  2©c)  (46) 

/?(q)  =  4i2(q)0o  (4c) 

As  readily  seen,  in  the  case  of  the  nonuniform  thin- 
walled  beam  theory,  rji,  Ac,  /?,  F^,  r„  and  a  are 


=  (66) 

^y(^.0  =  7*i(«,0-“o(^.0  (6c) 

where  Og  (z,t)  and  Oy  (z,t)  denote  the  rotations  about 
axes  z  and  y,  respectively,  7*,  and  7y,  denote  the 
transverse  sheeir  strains  in  the  planes  zz  and  yz,  re¬ 
spectively,  and  primes  denote  derivatives  with  respect 
to  z.  The  quantities  Ug  (z,  t) ,  Vg  (z,  /)  2uid  Ug  (z,  t)  de¬ 
noting  the  rigid-body  translations  along  z,  y  and  z 
axes,  respectively,  and  5,(z,t),  6yiz,t)  and  4>iz,t) 
denoting  the  rotations  about  z  and  y  axes  and  the 
twist  about  the  z-axis,  respectively,  represent  the  un¬ 
knowns  of  the  problem. 

Constitutive  Equations 

Consider  the  case  of  composite  thin-walled  beams 
consisting  of  a  finite  number  N  of  homogeneous  lay¬ 
ers.  It  is  assumed  that  the  material  of  each  con¬ 
stituent  layer  is  linearly  elastic  and  anisotropic  and 
that  the  bonding  between  the  layers  is  perfect.  The 
three-dimensional  constitutive  equations  for  a  gener¬ 
ally  orthotropic  elastic  material  can  be  expressed  in 
the  matrix  form 


■c-„  ■ 

■Qii 

Qi7 

Qi3 

0 

0 

Ois" 

O12 

Q22 

Q23 

0 

0 

Q26 
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Cnn 
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_0 

Oas 

^nn 
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Q44 

Qas 
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Izn 
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Qss 
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where  Qij  denote  the  transformed  elastic  coefficients 
associated  with  the  ibth  layer  in  the  global  coordinate 
system  of  the  structure  and  7pr  =  2epr ,  ^  ^ij 

denote  the  components  of  the  strain  tensor.  The  three- 
dimensional  dependence  in  Eqs.  (7)  can  be  reduced 
to  an  equivalent  one-dimensional  dependence  in  two 
steps.  The  first  step,  yielding  the  two-dimensional 
local  constitutive  equations,  consists  of  the  integration 
of  the  original  three-dimensional  form  through  the 
Iciminate  thickness,  while  the  second  step,  resulting  in 
the  one-dimensional  form,  consists  of  the  integration 
of  the  previous  form  of  constitutive  equations  along 
the  midline  contour  of  the  beam  cross-section. 

In  light  of  i)  and  ii),  and  considering  Eqs.  (5)  and 
(7),  the  local  consiiiuiivt  equations  expressed  in  terms 
of  the  strain  measures  can  be  shown  to  consist  of  the 
equations  for  the  strefi  resultants 


shear  strains  in  the  mid-surface  of  the  beam  induced  by 
transverse  shear  and  by  the  twist,  respectively,  and 
denotes  the  transverse  shear  strain  component.  The 
stress  resultants,  the  stress  couples  eind  the  strain  mea¬ 
sures,  Eqs.  (8)-(10),  exhibit  a  two-dimensional  spatial 
dependence,  the  dependence  being  on  s  amd  z.  In  the 
dynamical  problem,  they  also  depend  on  time.  More¬ 
over,  in  light  of  the  dependence  of  Ac  and  on  z,  the 
stiffness  coefficients  Kiz,  K23,  K43,  K53  depend  on  z, 
as  well. 

The  Boundary- Value  Problem 

The  boundary-value  problem,  consisting  of  the 
differential  equations  of  motion  and  the  boundary 
conditions,  can  be  derived  conveniently  by  means 
of  the  extended  Hamilton’s  principle,^  which  can  be 
stated  as  follows: 


'Nc/ 

'Kn 

Kx2 

K\3 

Ki,' 

K21 

K22 

K23 

K24. 

,0  n 


the  transverse  shear  stress  resultant 

and  the  equations  for  the  stress  couples 


'Lez' 

'K^i 

Ka2 

Ka3 
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(10) 


where  Kij  denote  the  modified  local  stiffness  coeffi¬ 
cients  (see  Appendix).  Consistent  with  Eqs.  (5),  the 
strain  components  entering  into  the  above  constitutive 
equations  are  e,,  =  -h  ncj,,  7,,  =  -f  and 
Tzn,  where 

f°^(s,z,t)  =u;o  -t-z(5,z)0y 

+  y{s,z)e',-F^{s,z)<f>"  (11a) 

c^,(s,z,t)=f?;^-5;g-a(s,z)<^"  (116) 

7sz{s,Z,t)  =[Uo  +  K  +^r]^  (lie) 

-j,,{s,z,t)=2^f>'  (lid) 

7zn{s,  Z,t)  =[Uo  +  “  K  +  (Hc) 

in  which  c®,  and  denote  the  axial  strain  compo¬ 
nents  associated  with  primary  and  secondary  warp¬ 
ing,  respectively,  7,^  and  7,,  represent  the  tangential 


r 


{6T  -  6V  +  6W)  dt  =z  0, 


6uo=zSvo~6wo  =  69x  —  66y  —  69  =  0  at  t  —  ti,t2  (12) 


where  T  is  the  kinetic  energy,  V  the  potential  energy, 
which  is  equal  to  the  strain  energy,  and  SW  the  virtual 
work  due  to  nonconservative  forces.  The  kinetic  energy 
has  the  form 


+  (vo  +  *(s.  +1  («J0  +  !/(s,  i)K 

-(-x(s,  z)9y  -  Fuis,  z)<j>'^  -f  n  ^^^y 


dx  A  ,  . 

-  a(s,z)(p 


dndsdz 


(13) 


On  the  other  hand,  the  strain  energy  can  be  shown  to 
have  the  expression 


-FUs,  z)4>"  +  n  -  a(s,  z)<6") 


U)o-t-i(s,z)0'  -t-y(s,r)0r 


(t) 


,  ,  ^  V  dy  ,  ,  „  ,  dz 

(Uo  +  0,)--(Uo  +  e.)-jj 


>  dndsdz  (14) 


Moreover,  the  virtual  work  of  nonconservative  forces 
can  be  written  as 


SW  zz  /  (pj,6t;o -I- m,6(6)dz 

Jo 


(15) 
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where  and  m,  denote  the  lift  force  per  unit  length 
and  aerodynamic  twist  moment  (positive  nose  up) 
about  the  elastic  axis. 

Carrying  out  integrations  with  respect  to  n,  s  and 
t,  we  can  write 


6Tdt  =  - 


(IiSuo  + 


+  {h  -  I9)  s<t>  +  hSdy  +  IrSe^)  dz  -  h<t> 


dt  (16) 


and  transverse  shear.  It  is  the  principal  goal  of  struc¬ 
tural  t^oring  to  select  the  appropriate  fiber  orienta¬ 
tion  so  as  to  produce  desired  elastic  couphngs  between 
certain  modes.  For  the  problem  at  hand,  the  induced 
elastic  couplings  must  play  a  decisive  role  m  enhanc¬ 
ing  the  free  vibration  and  aeroelastic  response  char¬ 
acteristics  of  the  wing  structure.  In  this  sense,  the 
bending-twist  cross-coupling  is  the  most  significant  in 
the  design  of  aircraft  wings.  The  above  criteria  for 
selecting  fiber  orientation,  together  with  ease  of  im¬ 
plementation  in  design  and  manufacturing,  result  in  a 
ply-angle  distribution  governed  by 


where  7,  denote  inertia  terms  (see  Appendix),  func¬ 
tions  of  z  and  t.  Upon  integration  with  respect  to  n 
and  s  in  Eq.  (14),  we  can  write  the  variation  in  the 
potential  energy  « 

SV  =  -  J  [TiSwo  +  {My  -  Qt)  ^6y 

+  (Af;  -  Q„)  -b  (5" -b  MO  5^  +  q;«uo 
-l-Qy(5vo]  dz  +  [TiSwo  -b  MyS6y  +  Mx66x  -  Bu64> 


+  (Bi  +  Mz)64>  +  Qzbuo  +  Qy6vo] 

0 

(17) 

where 

Tz{z,t)  =  j  Nzzds 

(18a) 

Qr{z,t)  =  j{N.z^  +  N,n-£jds 

(185) 

QA^,t)  =  f{N,zf^-Nzn"^)ds 

(18c) 

Mx(z,0  =  ^  (y{s,z)Nzz  - 

(18d) 

My{z,t)  =  ^  (^x{s,z)Nzz+Lzz^'^ 

(18c) 

Mz{z,t)  =  2  N,zip{z)ds 

(18/) 

Buiz,t)  =  [Fu{s,z)N„  +  a{s,z)L,z]ds 

(18*/) 

In  Eqs.  (18),  Tz,  Qz  and  Qy  denote  the  axial  force  and 
shear  forces  in  the  z-  and  y-directions,  M*,  My  and 

Mz  denote  the  moments  about  the  z-,  y-  and  r-axes, 
respectively,  and  is  the  bimoment  quantity.  Sub- 
stition  of  Eqs.  (15)-(17)  into  Eq.  (12),  followed  by  the 
usual  steps,  results  in  the  boundary-value  problem  for 
the  most  general  case  of  anisotropy.  The  boundary- 
value  problem  consists  of  six  differential  equations  of 
motion  with  variable  coefficients  for  the  displacements 
UQ,VQ,Wo,Bx,ey  together  with  the  correspond¬ 

ing  boundary  conditions.  Such  a  set  exhibits  complete 
coupling  between  the  various  modes,  i.e.,  primary  and 
secondary  warping,  vertical  and  lateral  bending,  twist 


B{y)  =  -B{-y)- 

The  ply-configuration  is  shown  in  Fig.  3.  According 
to  terminology  adopted  in  Refs.  8  and  9,  structures 
displaying  this  ply-angle  distribution  are  referred  to 
as  circumfenntially  asymmetric  stiffness  configuration 
and  symmetric  configuration,  respectively. 

As  a  result  of  the  ply-angle  distribution,  Eq.  (19), 
Hamilton’s  principle,  Eq.  (12),  yields  two  independent 
boundary-value  problems,  referred  to  here  as  Problem 
A  and  Problem  B.  Problem  A  is  of  eighth  order  and 
involves  the  twist  4>,  the  vertical  bending  vo  and  the 
vertical  transverse  shear  Bx-  On  the  other  hand, 
Problem  B  is  of  sixth  order  and  involves  the  extension 
luo,  the  lateral  bending  uq  and  the  lateral  transverse 
shear  By.  For  cantilevered  beams,  boundary- value 
Problem  A  is  governed  by  the  differential  equations 
of  motion 

-  (066^+  {^M'  +  (“77^')'  + 

=  ^^(64  -f  65)(ji  -  ((^10  +  he)  (20“) 

[055  (t^'o  +  «r)]'  +  (a56<^”)'  +  BdPy  =  (206) 

{asaB'x)'  +  {a374>')'  -  a55(i'o  +  ^z)  -  £56^ 

=  6v  (64  +  614)^*  (20c) 

to  be  satisfied  over  0  <  z  <  Moreover,  at  z  =  0 
the  solution  of  Eqs.  (20)  must  satisfy  the  boundary 
conditions 

^  =  0,  vo  =  0,  ffr  =  0,  ^  =  0  (21a -d) 

and  at  z  =  L  it  must  satisfy  the  boundary  conditions 

-  (a66<6")'  +  073^1  +  =  -6v(6io-b6i8)^'  (22a) 

+  Bz)  +  °56<6"  =  0  (22^) 

033^*  +  “37*6'  =  0,  066*6"  =  0  (22c,  d) 

Note  that  the  underlined  terms  in  Eqs.  (20)-(22)  are 
due  to  the  w2tfping  restraint  effect.  Moreover,  to  iden¬ 
tify  the  two  problems  to  be  studied,  namely,  the  free  vi¬ 
bration  and  the  aeroelastic  divergence  instability,  two 
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tracing  quantities,  6^  and  6d,  have  been  introduced  in 
Eqs.  (20)-(22).  For  Problem  A,  =  1  and  Sj  =  0, 
whereas  for  Problem  B,  =  1  and  =  0- 

The  boundary-value  Problem  B  is  governed  by  the 
differential  equations 

(a4iU'o)^  +  [a44(Wo  +  (23a) 

(aiiw'g)'  +  [oi4(tio  -I-  Sy)]  =  biiiio  (236) 

(a22^i)^  -  a4iiUo  -  a44(u'<,  4-  6y)  =  (65  -f  6i5)5y  (23c) 

and  the  boundary  conditions 

uo  =  0,  luo  =  0,  6y=0  (24a-c) 

to  be  satisfied  at  z  =  0,  as  well  as  the  boundary 
conditions 

+  “44(^0  -I-  ^y)  =  0 

aiiiWo  4- ai4(u(, -I- tfy)  =  0  (25) 

a220y  =  0 

to  be  satisfied  at  2  =  i. 

In  the  following,  the  study  of  the  boundary-value 
problem  A  governing  the  vibrational  and  aeroelastic 
response  behavior  of  wing  structures  will  be  carried 
out.  It  should  be  mentioned  here  that  the  stiffness 
terms  037  =  073  and  use  =  aes  appearing  in  Eqs.  (20) 
and  (22)  are  responsible  for  the  coupling  between 
bending  and  twist,  with  the  effect  of  055  expected 
to  be  weaker  than  the  effect  of  037.  The  boundary- 
value  problem  B  involves  only  a  single  coupUng  term, 
014  =  a4i ,  coupling  the  lateral  transverse  shear  and  the 
cixial  extension  motions.  The  rigidity  coefficients  o,; 
and  the  inertia  coefficients  6.-  appearing  in  Eqs.  (20), 
(22),  (23)  and  (25)  are  displayed  in  the  Appendix. 
With  the  exception  of  an  and  61,  all  the  other  rigidities 
Qij  and  inertia  coefficients  6,-  depend  on  2. 

Structural  Tailoring  for  Improved  Vibration 
aind  Static  Aeroelastic  Instability  Characteris¬ 
tics 

Static  aeroelastic  instability  is  an  important  fac¬ 
tor  in  the  design  of  modern  aircraft.  The  analysis  per¬ 
formed  here  addresses  the  problem  of  designing  the 
wing  so  as  to  take  advantage  of  structural  coupUngs 
from  a  static  aeroelastic  viewpoint.  This  is  done  by 
using  the  unique  directional  properties  of  advanced 
composite  materials.  The  same  importance  should  be 
afforded  to  the  vibrational  characteristics,  which  are 
basic  to  determining  the  dynamic  response  and  flutter 
instability  and  to  actively  controlling  the  wing  struc¬ 
ture. 

In  the  case  of  free  vibration,  the  terms  associ¬ 
ated  with  the  external  loadings  are  omitted,  6^  =  0 


and  6v  =  1,  whereas  for  static  aeroelastic  problems, 
the  only  loading  terms  to  be  retained  are  the  ones 
associated  with  the  aerodynamic  hft  py  and  the  tor¬ 
sional  aerodynamic  moment  m, ,  amounting  to  6^  =  1 
and  6v  =  0.  Using  strip-theory  aerodynamics,  we  can 
write^° 

Py  (2)  =  qnc{z)ao  [<6o  4-  ^  -  Vg  tan  A]  -  NW/2L  (25o) 
m,{z)  =  qnc{z)age{z)  [4)e  +  <i>-  I'o  tan  A] 

+  qnciz)^CMAC  -  NWd/2L  (256) 

where  denotes  the  dynamic  pressure  nor¬ 

mal  to  the  leading  edge  of  the  swept  wing,  c(z)  the 
chord  of  the  wing,  Oo  the  “corrected  lift”  curve  slope  co¬ 
efficient,  A  the  angle  of  sweep  (considered  positive  for 
swept-back),  e(2)  the  of&et  between  the  aerodyn2unic 
and  reference  axis,  <j>o  the  rigid  angle  of  attack  (mea¬ 
sured  in  planes  normal  to  the  leading  edge),  Cmac  the 
wing  section  pitching  moment  coefficient  (whose  influ¬ 
ence,  as  usual,  is  disregarded),  WI2L  the  wing  weight 
per  unit  length  and  N  the  load  factor  normaJ  to  the 
wing  surface. 

Due  to  the  relatively  high  order  of  the  problem 
and  the  fact  that  the  problem  is  characterized  by  vari¬ 
able  coefficients,  the  solution  tends  to  be  very  complex. 
Towards  the  goal  of  solving  the  two  eigenvalue  prob¬ 
lems  a  powerful  computational  methodology  based  on 
the  extended  Galerkin  method  was  devised.  For  some 
special  cases,  the  solution  accuracy  was  checked  by 
comparing  the  results  with  the  exact  ones  obtained  via 
the  Laplace  transform  method,  devised  by  the  sernie 
authors,  and  the  agreement  was  found  to  be  excellent. 
As  a  generEil  remark,  we  observe  from  Eqs.  (25)  that 
for  A  <  0,  i.e.,  for  swept-forward  wings,  the  aeroe¬ 
lastic  bending-twist  coupling  results  in  an  increase  in 
Py(2)  and  m,(z),  which  in  turn  reduces  dramatically 
the  divergence  speed. 

The  study  of  the  static  aeroelastic  critical  case 
implies  determination  of  divergence  instabihty  condi¬ 
tions.  The  determination  of  the  divergence  speed  leads 
to  the  solution  of  an  eigenvalue  problem,  where  the  di¬ 
vergence  speed  plays  the  role  of  an  eigenvalue.  Struc¬ 
tural  tailoring  applied  to  the  vibration  of  wing  struc¬ 
tures  must  result  in  an  increase  in  the  eigenfrequencies 
without  weight  penalties.  The  determination  of  nat¬ 
ural  frequencies  requires  the  solution  of  an  eigenvzdue 
problem  as  well.^ 

Uncoupled  Cases  for  Divergence  Instability 

In  two  uncoupled  cases,  closed-form  solutions  for 
the  divergence  instabihty  can  be  obtained.  They 
correspond  to  a)  pure  bending  divergence  of  swept 
wings  infinitely  rigid  in  transverse  shear  and  b)  pure 
torsionzd  divergence.  In  the  former  case,  ehminating 
056(6"  from  the  Eqs.  (20b,c),  assuming  a  very  large 
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torsional  stiffness,  letting  Bx  — ►  —v'o  and  =  0,  id  =  1 
and  implementing  a  Rayleigh-quotient  procedure  to 
the  resulting  equation,  the  expression  for  the  obtained 
divergence  speed  can  be  shown  to  be 

-2  /  033(2)  {v'of  dz 

(9n)z,  = - - 7-  (26) 

Oo  tanA  /  c{z)  (vg)'  dz 

This  equation  revecils  that  only  a  swept-forward  wing 
(A  — >  -A)  can  exhibit  divergence  instability  in  pure 
bending.  This  result  constitutes  an  extension  to  the 
case  of  wings  modeled  as  thin-walled  beams  of  that 
obtained  in  Ref.  3  for  solid  beams. 

In  the  pure  torsional  case,  assuming  that  the 
bending  stiffness  is  large  and  implementing  a  Rayleigh 
quotient  procedure  in  conjunction  with  Eq.  (20a),  one 
obtains 


(9n)D  =qD  = 

dz+  {  (077(2))  (4>'f  dz 

h - - - h -  (27) 

Go  I  c(z)e{z)<fy^dz 
Jo 

in  which  the  underscored  term  is  connected  with  the 
warping  inhibition.  As  in  the  case  of  solid  beams,® 
Eq.  (27)  reveals  that  pure  torsional  divergence  can 
occur  for  straight  wings  only. 

From  Eq.  (26),  we  conclude  that,  to  increase 
as  much  as  possible  {qn)o,  the  wing  structure  in 
pure  bending  should  be  designed  so  that  K14  —>■  0, 
a  condition  resulting  in  an  increase  of  the  overzdl 
bending  stiffness  033. 


Closed-Form  Solution  for  the  Divergence  of 
Swept  Forwrurd  Wings 


Having  the  conditions  corresponding  to  decoupled 
divergence  in  bending  and  torsion,  Eqs.  (26)  and  (27), 
respectively,  a  closed-form  solution  for  the  coupled 
divergence  of  swept  forward  wings  can  be  established. 
Although  approximate,  for  the  case  of  solid  beams,^ 
the  solution  was  shown  to  provide  results  in  excellent 
agreement  with  the  exact  ones. 

The  solution  for  the  divergence  of  swept  forward 
wings  is  based  on  a  linear  algebraic  relationships  of  the 
two  decoupled  expressions  of  the  divergence  instability. 
The  linear  relationship  yields  the  expression  for  the 
divergence  in  the  form 


{qn)D 


1  U 
OgCjlL^  S 


(28) 


where 


U  = 


i: 

Aee  066(2) 
077(2) 

(ff) 

2' 

1 

df] 

j 

CE<f>^dT) 

0 

(29a) 

and 

divl) 

t/ItanA  /  C-^dr) 

S  =  - ^ ,  (296) 

In  Eqs.  (29),  Aee  =  066/(066)^)  -^77  =  ar7/(a77)fi 
and  j433  =  033/  (033)^  are  normalized  stiffness  quan¬ 
tities  while  C  =  c(z)/cfl  and  E  =  e{z)/eii  denote  the 
normalized  chord  emd  distance  between  aerodynamic 
line  and  reference  axis,  respectively.  It  should  be  men¬ 
tioned  that  in  Eqs.  (28),  (29)  the  warping  restraint 
effect  was  included  but  the  effect  of  transverse  sheu 
flexibility  was  not.  Equation  (28)  represents  the  coun¬ 
terpart  for  thin- walled  beams  of  the  expression  for  the 
divergence  obt2dned  in  Ref.  3  for  solid  beams. 

Niunerical  Results  and  Discussions 

The  wing  structure,  whose  geometrical  chairacter- 
istics  are  depicted  in  Fig.  1,  is  assumed  to  be  made 
of  graphite-epoxy  materizd  with  the  following  elastic 
properties: 

El  =  30  X  10®  psi,  Ez  =  E3  =  0.75  x  10®  psi 
Gi3  =  G23  =  0.37  X  10®  psi,  Gi2  =  0.45  x  10®  psi 

=  f^23  =  A^i3  =  0.25,  p  =  14.3  X  10~®  lb  sec^/in^ 

The  eigenvalue  problem  was  solved  for  a  variety  of 
cases  and  the  results  are  presented  in  the  following. 
Figures  3  and  4  show  the  first  three  eigenfrequencies 
versus  the  ply  angle  for  various  values  of  the  taper 
ratio.  Figures  5  emd  6  display  the  first  four  uncoupled 
eigenfrequencies  versus  the  taper  ratio  for  6  =  0°  2uid 
6  =  90'  and  Fig.  7  presents  the  first  four  coupled 
eigenfrequencies  for  6  =  45' .  Figures  8  and  9  show  the 
distribution  of  the  normalized  vertical  displacement 
2ind  twist  in  the  lowest  three  modes  for  6  =  45'  and 
for  different  values  of  taper  ratio.  Figure  10  displays 
the  divergence  speed  versus  the  taper  ratio  for  a  swept 
back  wing  with  A  =  5'  and  with  ply  angles  0  =  0'  and 
90'.  Finally,  Figs.  11  and  12  display  the  distribution 
along  the  wing  span  of  the  normalized  bending  and 
torsional  stiffness.  The  figures  reveal  that  the  linear 
variation  along  the  span  of  the  chord  and  height  of  the 
cross-section  profile  induces  a  pwabolic  distribution 
for  these  stiffness  quantities. 
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In  all  the  plots  presented,  the  significant  effects  of 
the  taper  ratio  on  the  vibration  and  static  aeroelastic 
instability  characteristics  of  aircraft  wings  have  been 
demonstrated.  As  shown,  an  increase  in  the  taper 
ratio  results  in  a  decay  in  the  vibrational  and  static 
aeroelastic  properties  of  the  wing  structure.  This 
trend  can  be  reversed  by  implementation  of  tailoring 
techniques.  This  amounts  to  proper  selection  of  the 
ply  angle  so  as  to  reduce  negative  effects  of  wing  taper. 
Hence,  tailoring  permits  exploitation  of  advantages 
offered  by  wing  taper. 
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APPENDIX 

The  modified  local  stiffness  coefficients  are  given 
by 


Kn  =  Ai,  -  4^,  Ki2  =  A26  -  =  K21 

An  All 

KM  =  ^Kn^>  Kn  =  B22  - 

K22  =  Aas  -  4^,  K23{z)  =  2K22^ 

All 


(A.l) 


^  n  AieBia  _ 

■ft  24  —  -026 - -  —  ^42 

All 

KM  =  ^K24^,  =  D22-^^ 

Tj-  _  n  -SisAia  _  D  -BieAis 

A51  =  1326 - - .  A52  =  £>66 - - 

All  All 

KM  =  2K62^,  Kn  =  D26  - 

P(Z)  All 

where  A^,  Bij  and  Dij  denote  local  stretching, 
stretching-bending  coupling  and  bending  rigidity 
quantities,  respectively.  The  dependence  on  the  z- 
coordinate  of  various  coefficients  components  Kij  is 
indicated  explicitly. 

The  inertia  terms  are  as  follows: 
h{z)=^  [fio  -  y(«)^]  mods 
h  =j>  [vo  +  mods 

Isiz)  [ino  +  x{z)6y  +  y(z)Pr  -  Fuis,  z)^'  mods, 

U{z)=  +  y^{z))4>  -  y('j)«o  +  a:(2:)i;o]rnods 

^5(^)  [z(^)u)o  +  x^{z)ey 

+  x{z)y{z)6x  -  x{z)Fu{s,  z)^']mods 

/(A. 2) 

[y(z)uio  +  y^{z)6x 

+  x{z)y{z)ey  -  y{z)Fu{8,  z)^']mod« 
f  \dxdx  "  dx  dy  ■■  dx  ,  .  ••,1 

hi!)  =  f[-Fu[s,z)wo  -  x{z)Fu{i,z)9, 

-  y{z)Fu{s,  z)6x  -I-  Fuis,  z)^(i)']mods 
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+a^(5,z)(^']  m2ds 


where 


(mo,m2)  =  Y^/  P(k)(i,f^^)dn  (^.3) 

fc_i  JA(k_,) 

denote  mass  terms. 

The  rigidity  coefficients  have  the  form 
On  =  y  Knds,  012(2)  =  ^  |^a:(2)A’ii  +  ^14^ 
ai3(^)  =  ^ 

/dx  f  dy 

Ki2-^ds,  015(2)  =  ®  Ki2—ds, 

016(2)  =  ~  ^  [■f!^ii'f’w(s.  2)  +  Ki4ais,  z)]  ds 

017(2)  =  ^  Ki3{z)ds 

022(2)=  [/£:iix^(2)  +  2x(2)/i:i4;;57  +  /f44;;n;r 


ds. 


023(2) 


=/[ 


'ds 

Knx{z)yiz)  -  x{z)Ki4^ 


'dsds 


,  dy  dxdy 

+y(2)/^14^-/^44^jj 


ds 


024(2)  =  y  ( 
025(2)  =  J  ^ 


x(2)Xi2^  +  i^24jf)<i5 


xiz)K,2^  +  K24Y,^ys 


026(2)  =  -  ^  [x{z)KuF^{s,  z)  +  x{z)Ki4a{s,z) 


+Fuis,z)Ki4-J^  +  K44a{8,z)-^ 


ds 


027(2)  =  y  |x(2)ii:i3(2)  +  K43{z)^  ds 
033(2)  =  y  j^i^iiy^2)  -  2y{z)Ki4^ 
034(2)  =  ^ 

035(2)  =  y 

036(2)  =  -  ^[yiz)KnF^+yKi4a{s,z) 


dx  dx  dx 

I  «^y  r.  dxdy 

y(2)i^i2--/(:24jj^ 


ds 

ds 


(AA) 


-Fu{s,z)Ki4^  -  K4ia{s,z)^ 


ds 


037(2) 


=/[ 


yiz)Ki3{z)  -  ■f^43(2)^ 


ds 


045(2)  =  -  y  5 1. 

A-r'i 

ds 


F^is,z)K2i^  +  K2Ms,z)^ 


ds 


047(2)  =  ^ 
055(2)  = 


dydy  dx  dx 

^^^TsTs-^^^TsTs 


F„(s,  z)K2i  ^  +  K24a{s,  2)^ 


ds 


056(2)  =  -^  J. 

057(2)  =  y 

066(2)  =  ^  [F'iiF*(s,  2)  +  2Ki4Fu{s,  z)a{s,  2) 
+A’44a*(s,  2)]  ds 

057(2)  =  -  y  [^13(2)Fu,(5,  2)  +  K43{z)a(s,  2)]  ds 


077(2) 


K23{z) 


ds 


and  note  that  the  rigidity  coefficients  are  synometric, 

o«;  ~  Oji(»,  j  —  1, 2, ...  7). 

Finally,  the  inertia  coefficients  are  given  by 


[61,  64(2),  65(2),  610(2)] 

=  ^mo[l,  y*(2),  1^(2),  F^i8,z)]ds 


[614(2),  615(2),  618(2)] 


ds 


(^.5) 


X 


i  =80in: 

en  siOin;  i 

bn  »  2m: 

h  at  0.4iii  { 

Fig.  la  Top  view  of  the  tapered  wing 
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